Rio Grande
Brian
1/30/2023
Disclaimer
This work is very preliminary as I get back into the coding swing of things. Data wrangling and figure generation will be done via R, but the rest of the project will be done using good ol’ microsoft products. This is just an entry point into data crunching and should by no means be considered a final product.
Also, I’m not great at this but whatever. I could automate this, but I’ll figure that out shortly!
Need to update figures & analyze sites for seasonal trends
Methodology
SNOTEL data was provided by the NRCS. Data was cleaned by removing outliers that are likely implausible; any year with more than 15 observations missing was removed. Temperatures were adjusted using the Morrisey method for stations identified by Ma et al (2019) due to SNOTEL temperature sensor changes, with the adjustment applied to pre-sensor change data. Daily mean observations were detrended to determine whether values were increasing or decreasing from the entire time series trend. Daily mean temperatures were first averaged by water year, with all water year means then averaged by day of water year. The mean temperature by day for the period of record was averaged. To find the standard deviation, the daily mean temperatures by water year was subtracted from the averaged mean temperature by day for the period of record. All water year means averaged by day of water year were subtracted from the temperature mean. The resulting values were then added together to find the “residual” of the daily mean temperatures by water year. The standard deviation was then computed from those residuals, with trends analyzed by Mann‐Kendall significance test and Theil‐Sen’s rate of change. Significant trends are identified with p-values of less than 0.10.
Morrisey Method
The Morrisey Method is taken from Ma, Fassnacht and Kampf..
In R script: T(adjusted) = 5.3x10(-7)xT(old)4+3.72x10(-5)xT(old)3-2.16x10(-3)xT(old)2-7.32x10^(-2)xT(old)+1.37
In the Ma et al. spreadsheet, H1 is Morrisey, H2 is Oiler
San Juan Area SNOTEL sites:
Beartown 327 Original 4/18/2005
Cascade #2 389 Morrisey 6/18/2004
Cumbres Trestle 431 Oyler -> Morrisey 6/8/2005
Idarado 538 Morrisey 7/12/2004
Lone Cone 589 Oyler -> Morrisey 6/22/2005
Middle Creek 624 Morrisey 6/26/2006
Mineral Creek 629 Oyler ?? -> Morrisey 6/22/2004
Molas Lake 632 NSCE-Morrisey, Bias-Original 10/2/2003
Red Mountain Pass 713 Morrisey 8/18/2004
Scotch Creek 739 Morrisey 6/15/2004
Slumgullion 762 Morrisey 6/26/2006
Spud Mountain 780 Morrisey 6/28/2004
Stump Lakes 797 Morrisey 7/22/2005
Upper Rio Grande 839 Oyler -> Morrisey 5/26/2004 *Why is this in red?
Upper San Juan 840 Morrisey 4/7/2004
Vallecito 843 Morrisey 7/22/2005
Wolf Creek Summit 874 Morrisey 7/12/2004
Data from thesis research:
SNOTEL_san_juan_Area <- snotel_download(site_id = c(327, 389, 431, 538, 589, 624, 629, 632, 713, 739, 762, 780, 797, 839, 840, 843, 874), path = tempdir('../data'), internal = TRUE)
write.csv(SNOTEL_san_juan_Area,"C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", row.names = FALSE) #write in the raw data
Beartown 327
Original
SNOTEL_san_juan_Area <- read.csv("C:/Users/13074/Documents/ESS580/thesis_project/San_Juan_area/data_raw/snotel_san_juans.csv", header = TRUE)
snotel_327 <- SNOTEL_san_juan_Area %>%
filter(site_id == "327")#str(snotel_327) # check the date, usually a character.
snotel_327$Date <- as.Date(snotel_327$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_327_clean <- snotel_327 %>% # filter for the timeframe
filter(Date >= "1983-08-10" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_327_clean <- snotel_327_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_327_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_327_clean <- snotel_327_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40)ggplot(snotel_327_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_327_cull_count <- snotel_327_clean %>%
filter(temperature_min < -40) %>%
count(waterYear)
snotel_327_cull_count## # A tibble: 0 x 2
## # Groups: waterYear [0]
## # ... with 2 variables: waterYear <dbl>, n <int># filtering for too few observations in a year
snotel_327_cull_count_days <- snotel_327_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_327_cull_count_days## # A tibble: 3 x 2
## # Groups: waterYear [3]
## waterYear n
## <dbl> <int>
## 1 1983 52
## 2 1994 337
## 3 2003 3461983, 1994, 2003 need to be culled.
snotel_327_clean_culled <- snotel_327_clean %>%
filter(waterYear != "1983" & waterYear != "1994" & waterYear != "2003") #%>%
#filter(temperature_mean > -49)Beartown (327) NOT ADJUSTED.
2005-04-18 installed ext range sensor. Original
snotel_327_adjusted <- snotel_327_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-04-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
327 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_327 <- snotel_327_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_327 <- yearly_wy_aver_327 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_327 <- daily_wy_aver_327 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_327$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_327 <-daily_wy_aver_327 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_327$date_temp <- signif(daily_wy_aver2_327$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_327, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
327 SD
standard_dev_327 <- daily_wy_aver_327 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_327 <- standard_dev_327 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_327 <- standard_dev_all_327 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_327 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 8.215357 |
| 1985 | 8.122497 |
| 1986 | 7.129599 |
| 1987 | 7.619777 |
| 1988 | 8.257318 |
| 1989 | 8.044976 |
| 1990 | 7.675963 |
| 1991 | 7.684970 |
| 1992 | 7.079517 |
| 1993 | 7.642300 |
| 1995 | 7.304999 |
| 1996 | 7.681058 |
| 1997 | 7.711950 |
| 1998 | 8.105700 |
| 1999 | 6.825322 |
| 2000 | 7.526413 |
| 2001 | 7.912097 |
| 2002 | 8.099855 |
| 2004 | 7.622635 |
| 2005 | 7.225545 |
| 2006 | 7.107924 |
| 2007 | 7.426444 |
| 2008 | 7.802192 |
| 2009 | 6.893765 |
| 2010 | 7.820347 |
| 2011 | 7.632307 |
| 2012 | 7.466342 |
| 2013 | 8.000100 |
| 2014 | 7.260505 |
| 2015 | 6.671562 |
| 2016 | 7.474460 |
| 2017 | 7.072348 |
| 2018 | 6.933787 |
| 2019 | 7.658054 |
| 2020 | 7.769482 |
| 2021 | 7.841403 |
| 2022 | 7.318492 |
ggplot(standard_dev_all_327, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 327 average temperatures for water years 2005-2021
MK & SS for 327 (non-corrected)
sd_mk_327 <- mk.test(standard_dev_all_327$sd_2)
print(sd_mk_327)##
## Mann-Kendall trend test
##
## data: standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -140.0000000 5846.0000000 -0.2102102sd_sens_327 <- sens.slope(standard_dev_all_327$sd_2)
print(sd_sens_327)##
## Sen's slope
##
## data: standard_dev_all_327$sd_2
## z = -1.818, n = 37, p-value = 0.06907
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.025461316 0.001275883
## sample estimates:
## Sen's slope
## -0.01223021Cascade #2 389
Morrisey 6/18/2004
snotel_389 <- SNOTEL_san_juan_Area %>%
filter(site_id == "389")#str(snotel_389) # check the date, usually a character.
snotel_389$Date <- as.Date(snotel_389$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_389_clean <- snotel_389 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_389_clean <- snotel_389_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_389_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_389_clean <- snotel_389_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40) %>%
filter(temperature_mean > -40)ggplot(snotel_389_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
temp_389_xts <- xts(snotel_389_clean$temperature_mean, order.by = snotel_389_clean$Date)
dygraph(temp_389_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_389_cull_count <- snotel_389_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_389_cull_count## # A tibble: 40 x 2
## # Groups: waterYear [40]
## waterYear n
## <dbl> <int>
## 1 1983 325
## 2 1984 350
## 3 1985 350
## 4 1986 365
## 5 1987 358
## 6 1988 366
## 7 1989 330
## 8 1990 348
## 9 1991 364
## 10 1992 366
## # ... with 30 more rows# filtering for too few observations in a year
snotel_389_cull_count_days <- snotel_389_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_389_cull_count_days## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1983 325
## 2 1989 330
## 3 1990 348
## 4 1994 310
## 5 2006 330
## 6 2010 334
## 7 2014 314
## 8 2015 185
## 9 2016 3341983, 1989, 1990, 1994, 2006, 2010, 2014, 2015, 2016 need to be culled.
snotel_389_clean_culled <- snotel_389_clean %>%
filter(waterYear != "1983" & waterYear != "1990" & waterYear != "1994" & waterYear != "2006" & waterYear != "2010" & waterYear != "2014" & waterYear != "2015" & waterYear != "2016") #%>%
#filter(temperature_mean > -49)Cascade #2 (389) NOT ADJUSTED.
Morrisey 6/18/2004
snotel_389_adjusted <- snotel_389_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-06-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
389 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_389 <- snotel_389_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_389 <- yearly_wy_aver_389 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_389 <- daily_wy_aver_389 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_389$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_389 <-daily_wy_aver_389 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_389$date_temp <- signif(daily_wy_aver2_389$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_389, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
389 SD
standard_dev_389 <- daily_wy_aver_389 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389 <- standard_dev_389 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_389 <- standard_dev_all_389 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_389 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 10.776321 |
| 1985 | 7.879611 |
| 1986 | 10.887499 |
| 1987 | 9.408888 |
| 1988 | 10.711609 |
| 1989 | 10.418326 |
| 1991 | 10.101153 |
| 1992 | 8.637142 |
| 1993 | 8.884440 |
| 1995 | 9.371322 |
| 1996 | 9.429341 |
| 1997 | 10.027003 |
| 1998 | 9.698491 |
| 1999 | 9.029327 |
| 2000 | 9.288091 |
| 2001 | 10.630687 |
| 2002 | 10.346726 |
| 2003 | 9.727202 |
| 2004 | 9.282545 |
| 2005 | 8.472227 |
| 2007 | 9.686674 |
| 2008 | 9.541070 |
| 2009 | 8.659960 |
| 2011 | 9.334867 |
| 2012 | 9.418612 |
| 2013 | 9.751020 |
| 2017 | 9.144382 |
| 2018 | 8.911243 |
| 2019 | 9.427120 |
| 2020 | 9.880776 |
| 2021 | 9.863539 |
| 2022 | 9.441803 |
ggplot(standard_dev_all_389, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 389 average temperatures for water years 2005-2021
MK & SS for 389 (non-corrected)
sd_mk_389 <- mk.test(standard_dev_all_389$sd_2)
print(sd_mk_389)##
## Mann-Kendall trend test
##
## data: standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -60.0000000 3802.6666667 -0.1209677sd_sens_389 <- sens.slope(standard_dev_all_389$sd_2)
print(sd_sens_389)##
## Sen's slope
##
## data: standard_dev_all_389$sd_2
## z = -0.95677, n = 32, p-value = 0.3387
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.04702986 0.01336703
## sample estimates:
## Sen's slope
## -0.01646233Corrected
389 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_389_ad <- snotel_389_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_389_ad <- yearly_wy_aver_389_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_389_ad <- daily_wy_aver_389_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_389_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_389_ad <-daily_wy_aver_389_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_389_ad$date_temp_ad <- signif(daily_wy_aver2_389_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_389_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
389 SS (corrected)
standard_dev_389_ad <- daily_wy_aver_389_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_389_ad <- standard_dev_389_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_389_ad <- standard_dev_all_389_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_389_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 10.061639 |
| 1985 | 7.175032 |
| 1986 | 9.930924 |
| 1987 | 8.682951 |
| 1988 | 9.868631 |
| 1989 | 9.733860 |
| 1991 | 9.472969 |
| 1992 | 8.002025 |
| 1993 | 8.283300 |
| 1995 | 8.657379 |
| 1996 | 8.697682 |
| 1997 | 9.311693 |
| 1998 | 8.929771 |
| 1999 | 8.328247 |
| 2000 | 8.510889 |
| 2001 | 9.808164 |
| 2002 | 9.563945 |
| 2003 | 8.940614 |
| 2004 | 8.587303 |
| 2005 | 8.472294 |
| 2007 | 9.686174 |
| 2008 | 9.541322 |
| 2009 | 8.659292 |
| 2011 | 9.334237 |
| 2012 | 9.418205 |
| 2013 | 9.750812 |
| 2017 | 9.144257 |
| 2018 | 8.911176 |
| 2019 | 9.426664 |
| 2020 | 9.880403 |
| 2021 | 9.863015 |
| 2022 | 9.441185 |
ggplot(standard_dev_all_389_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 389 average temperatures for water years 1986-2021
MK & SS 389 (corrected)
sd_mk_389_ad <- mk.test(standard_dev_all_389_ad$sd_2)
print(sd_mk_389_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 64.0000000 3802.6666667 0.1290323sd_sens_389_ad <- sens.slope(standard_dev_all_389_ad$sd_2)
print(sd_sens_389_ad)##
## Sen's slope
##
## data: standard_dev_all_389_ad$sd_2
## z = 1.0216, n = 32, p-value = 0.307
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01140359 0.04896578
## sample estimates:
## Sen's slope
## 0.01785583Cumbres Trestle 431
Oyler -> Morrisey 6/8/2005
snotel_431 <- SNOTEL_san_juan_Area %>%
filter(site_id == "431")#str(snotel_431) # check the date, usually a character.
snotel_431$Date <- as.Date(snotel_431$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_431_clean <- snotel_431 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_431_clean <- snotel_431_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_431_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_431_clean <- snotel_431_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_min > -40) %>%
filter(temperature_min < 25)ggplot(snotel_431_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_431_cull_count <- snotel_431_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_431_cull_count## # A tibble: 36 x 2
## # Groups: waterYear [36]
## waterYear n
## <dbl> <int>
## 1 1987 1
## 2 1988 341
## 3 1989 364
## 4 1990 348
## 5 1991 341
## 6 1992 312
## 7 1993 356
## 8 1994 317
## 9 1995 364
## 10 1996 364
## # ... with 26 more rows# filtering for too few observations in a year
snotel_431_cull_count_days <- snotel_431_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_431_cull_count_days## # A tibble: 7 x 2
## # Groups: waterYear [7]
## waterYear n
## <dbl> <int>
## 1 1987 1
## 2 1988 341
## 3 1990 348
## 4 1991 341
## 5 1992 312
## 6 1994 317
## 7 2014 2881987, 1988, 1990, 1991, 1992, 1994, 2014 need to be culled.
snotel_431_clean_culled <- snotel_431_clean %>%
filter(waterYear != "1987" & waterYear != "1988" & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)ggplot(snotel_431_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)
dygraph(temp_431_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_431_clean_culled <- snotel_431_clean_culled %>%
filter(temperature_mean < 20)
temp_431_xts <- xts(snotel_431_clean_culled$temperature_mean, order.by = snotel_431_clean_culled$Date)
dygraph(temp_431_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Cumbres Trestle 431
Oyler -> Morrisey 6/8/2005
snotel_431_adjusted <- snotel_431_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-08", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
431 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_431 <- snotel_431_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_431 <- yearly_wy_aver_431 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_431 <- daily_wy_aver_431 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_431$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_431 <-daily_wy_aver_431 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_431$date_temp <- signif(daily_wy_aver2_431$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_431, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
431 SD
standard_dev_431 <- daily_wy_aver_431 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431 <- standard_dev_431 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_431 <- standard_dev_all_431 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_431 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1989 | 8.484011 |
| 1993 | 8.101376 |
| 1995 | 7.750446 |
| 1996 | 7.952021 |
| 1997 | 8.217323 |
| 1998 | 8.600343 |
| 1999 | 7.226028 |
| 2000 | 8.164156 |
| 2001 | 8.669911 |
| 2002 | 9.005161 |
| 2003 | 8.459072 |
| 2004 | 8.306872 |
| 2005 | 8.304679 |
| 2006 | 7.417238 |
| 2007 | 7.937180 |
| 2008 | 8.259566 |
| 2009 | 7.296965 |
| 2010 | 8.418225 |
| 2011 | 8.231369 |
| 2012 | 8.081428 |
| 2013 | 8.537331 |
| 2015 | 7.219164 |
| 2016 | 7.791453 |
| 2017 | 7.322790 |
| 2018 | 7.113402 |
| 2019 | 7.880447 |
| 2020 | 8.046531 |
| 2021 | 8.095523 |
| 2022 | 7.621146 |
ggplot(standard_dev_all_431, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 431 average temperatures for water years 2005-2021
MK & SS for 431 (non-corrected)
sd_mk_431 <- mk.test(standard_dev_all_431$sd_2)
print(sd_mk_431)##
## Mann-Kendall trend test
##
## data: standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -104.0000000 2842.0000000 -0.2561576sd_sens_431 <- sens.slope(standard_dev_all_431$sd_2)
print(sd_sens_431)##
## Sen's slope
##
## data: standard_dev_all_431$sd_2
## z = -1.9321, n = 29, p-value = 0.05335
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.041299369 0.001003281
## sample estimates:
## Sen's slope
## -0.01890302Corrected
431 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_431_ad <- snotel_431_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_431_ad <- yearly_wy_aver_431_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_431_ad <- daily_wy_aver_431_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_431_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_431_ad <-daily_wy_aver_431_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_431_ad$date_temp_ad <- signif(daily_wy_aver2_431_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_431_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
431 SS (corrected)
standard_dev_431_ad <- daily_wy_aver_431_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_431_ad <- standard_dev_431_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_431_ad <- standard_dev_all_431_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_431_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1989 | 7.914872 |
| 1993 | 7.518415 |
| 1995 | 7.173944 |
| 1996 | 7.365597 |
| 1997 | 7.634730 |
| 1998 | 7.970563 |
| 1999 | 6.682639 |
| 2000 | 7.537564 |
| 2001 | 8.039524 |
| 2002 | 8.344015 |
| 2003 | 7.813684 |
| 2004 | 7.737034 |
| 2005 | 7.599212 |
| 2006 | 7.418573 |
| 2007 | 7.938448 |
| 2008 | 8.260644 |
| 2009 | 7.298209 |
| 2010 | 8.419846 |
| 2011 | 8.232506 |
| 2012 | 8.082571 |
| 2013 | 8.538590 |
| 2015 | 7.220737 |
| 2016 | 7.793079 |
| 2017 | 7.324026 |
| 2018 | 7.114853 |
| 2019 | 7.881605 |
| 2020 | 8.047598 |
| 2021 | 8.096826 |
| 2022 | 7.622159 |
ggplot(standard_dev_all_431_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 431 average temperatures for water years 1986-2021
MK & SS 431 (corrected)
sd_mk_431_ad <- mk.test(standard_dev_all_431_ad$sd_2)
print(sd_mk_431_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 54.0000000 2842.0000000 0.1330049sd_sens_431_ad <- sens.slope(standard_dev_all_431_ad$sd_2)
print(sd_sens_431_ad)##
## Sen's slope
##
## data: standard_dev_all_431_ad$sd_2
## z = 0.99418, n = 29, p-value = 0.3201
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01334615 0.03134200
## sample estimates:
## Sen's slope
## 0.01029449Idarado 538
Morrisey 7/12/2004
snotel_538 <- SNOTEL_san_juan_Area %>%
filter(site_id == "538")#str(snotel_538) # check the date, usually a character.
snotel_538$Date <- as.Date(snotel_538$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_538_clean <- snotel_538 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_538_clean <- snotel_538_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_538_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_538_clean <- snotel_538_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_min < 25)ggplot(snotel_538_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_538_cull_count <- snotel_538_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_538_cull_count## # A tibble: 36 x 2
## # Groups: waterYear [36]
## waterYear n
## <dbl> <int>
## 1 1987 355
## 2 1988 366
## 3 1989 364
## 4 1990 365
## 5 1991 337
## 6 1992 358
## 7 1993 346
## 8 1994 351
## 9 1995 364
## 10 1996 365
## # ... with 26 more rows# filtering for too few observations in a year
snotel_538_cull_count_days <- snotel_538_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_538_cull_count_days## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 1991 337
## 2 1993 3461991 and 1993 need to be culled.
snotel_538_clean_culled <- snotel_538_clean %>%
filter(waterYear != "1991" & waterYear != "1993") # & waterYear != "1990" & waterYear != "1991" & waterYear != "1992" & waterYear != "1994" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)ggplot(snotel_538_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)
dygraph(temp_538_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_538_clean_culled <- snotel_538_clean_culled %>%
filter(temperature_mean < 20)
temp_538_xts <- xts(snotel_538_clean_culled$temperature_mean, order.by = snotel_538_clean_culled$Date)
dygraph(temp_538_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Idarado 538
Morrisey 7/12/2004
snotel_538_adjusted <- snotel_538_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-07-12", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
538 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_538 <- snotel_538_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_538 <- yearly_wy_aver_538 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_538 <- daily_wy_aver_538 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_538$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_538 <-daily_wy_aver_538 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_538$date_temp <- signif(daily_wy_aver2_538$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_538, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
538 SD
standard_dev_538 <- daily_wy_aver_538 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538 <- standard_dev_538 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_538 <- standard_dev_all_538 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_538 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.006364 |
| 1988 | 8.601864 |
| 1989 | 8.632529 |
| 1990 | 8.392986 |
| 1992 | 7.829013 |
| 1994 | 10.169815 |
| 1995 | 7.658097 |
| 1996 | 7.925011 |
| 1997 | 7.934668 |
| 1998 | 8.455869 |
| 1999 | 7.252626 |
| 2000 | 7.992580 |
| 2001 | 8.428201 |
| 2002 | 8.883209 |
| 2003 | 8.309823 |
| 2004 | 8.066471 |
| 2005 | 7.228205 |
| 2006 | 7.518896 |
| 2007 | 7.890984 |
| 2008 | 8.117130 |
| 2009 | 7.048527 |
| 2010 | 8.062733 |
| 2011 | 7.859864 |
| 2012 | 7.787906 |
| 2013 | 8.410412 |
| 2014 | 7.465408 |
| 2015 | 6.906228 |
| 2016 | 7.611126 |
| 2017 | 7.235848 |
| 2018 | 7.239780 |
| 2019 | 7.718304 |
| 2020 | 7.976082 |
| 2021 | 8.027478 |
| 2022 | 7.664387 |
ggplot(standard_dev_all_538, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 538 average temperatures for water years 2005-2021
MK & SS for 538 (non-corrected)
sd_mk_538 <- mk.test(standard_dev_all_538$sd_2)
print(sd_mk_538)##
## Mann-Kendall trend test
##
## data: standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -181.0000000 4550.3333333 -0.3226381sd_sens_538 <- sens.slope(standard_dev_all_538$sd_2)
print(sd_sens_538)##
## Sen's slope
##
## data: standard_dev_all_538$sd_2
## z = -2.6684, n = 34, p-value = 0.007621
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.042981556 -0.006896251
## sample estimates:
## Sen's slope
## -0.02384088Corrected
538 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_538_ad <- snotel_538_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_538_ad <- yearly_wy_aver_538_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_538_ad <- daily_wy_aver_538_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_538_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_538_ad <-daily_wy_aver_538_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_538_ad$date_temp_ad <- signif(daily_wy_aver2_538_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_538_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
538 SS (corrected)
standard_dev_538_ad <- daily_wy_aver_538_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_538_ad <- standard_dev_538_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_538_ad <- standard_dev_all_538_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_538_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.447584 |
| 1988 | 8.019241 |
| 1989 | 8.069873 |
| 1990 | 7.812557 |
| 1992 | 7.272544 |
| 1994 | 9.549437 |
| 1995 | 7.113641 |
| 1996 | 7.353811 |
| 1997 | 7.396458 |
| 1998 | 7.860736 |
| 1999 | 6.739299 |
| 2000 | 7.408528 |
| 2001 | 7.835100 |
| 2002 | 8.260161 |
| 2003 | 7.706375 |
| 2004 | 7.450714 |
| 2005 | 7.228187 |
| 2006 | 7.519259 |
| 2007 | 7.891254 |
| 2008 | 8.117349 |
| 2009 | 7.048722 |
| 2010 | 8.062844 |
| 2011 | 7.860215 |
| 2012 | 7.788573 |
| 2013 | 8.411054 |
| 2014 | 7.465688 |
| 2015 | 6.906947 |
| 2016 | 7.611742 |
| 2017 | 7.235817 |
| 2018 | 7.240024 |
| 2019 | 7.718445 |
| 2020 | 7.976612 |
| 2021 | 8.027693 |
| 2022 | 7.664991 |
ggplot(standard_dev_all_538_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 538 average temperatures for water years 1986-2021
MK & SS 538 (corrected)
sd_mk_538_ad <- mk.test(standard_dev_all_538_ad$sd_2)
print(sd_mk_538_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 1.100000e+01 4.550333e+03 1.960784e-02sd_sens_538_ad <- sens.slope(standard_dev_all_538_ad$sd_2)
print(sd_sens_538_ad)##
## Sen's slope
##
## data: standard_dev_all_538_ad$sd_2
## z = 0.14824, n = 34, p-value = 0.8821
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01508416 0.01703536
## sample estimates:
## Sen's slope
## 0.0009105465Lone Cone 589
Oyler -> Morrisey 6/22/2005
snotel_589 <- SNOTEL_san_juan_Area %>%
filter(site_id == "589")#str(snotel_589) # check the date, usually a character.
snotel_589$Date <- as.Date(snotel_589$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_589_clean <- snotel_589 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_589_clean <- snotel_589_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_589_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_589_clean <- snotel_589_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40)# %>%
#filter(temperature_min < 25)ggplot(snotel_589_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_589_cull_count <- snotel_589_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_589_cull_count## # A tibble: 37 x 2
## # Groups: waterYear [37]
## waterYear n
## <dbl> <int>
## 1 1986 1
## 2 1987 358
## 3 1988 366
## 4 1989 362
## 5 1990 365
## 6 1991 359
## 7 1992 366
## 8 1993 350
## 9 1994 337
## 10 1995 364
## # ... with 27 more rows# filtering for too few observations in a year
snotel_589_cull_count_days <- snotel_589_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_589_cull_count_days## # A tibble: 10 x 2
## # Groups: waterYear [10]
## waterYear n
## <dbl> <int>
## 1 1986 1
## 2 1994 337
## 3 1997 346
## 4 2006 317
## 5 2009 314
## 6 2010 332
## 7 2011 309
## 8 2012 322
## 9 2013 331
## 10 2014 3351986, 1994, 1997, 2006, 2009, 2010, 2011, 2012, 2013, 2014 need to be culled.
snotel_589_clean_culled <- snotel_589_clean %>%
filter(waterYear != "1986" & waterYear != "1994" & waterYear != "1997" & waterYear != "2006" & waterYear != "2009" & waterYear != "2010" & waterYear != "2011" & waterYear != "2012" & waterYear != "2013" & waterYear != "2014") #%>%
#filter(temperature_mean > -49)ggplot(snotel_589_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)
dygraph(temp_589_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_589_clean_culled <- snotel_589_clean_culled %>%
filter(temperature_mean < 20)
temp_589_xts <- xts(snotel_589_clean_culled$temperature_mean, order.by = snotel_589_clean_culled$Date)
dygraph(temp_589_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Lone Cone 589
Oyler -> Morrisey 6/22/2005
snotel_589_adjusted <- snotel_589_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-06-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
589 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_589 <- snotel_589_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_589 <- yearly_wy_aver_589 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_589 <- daily_wy_aver_589 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_589$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_589 <-daily_wy_aver_589 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_589$date_temp <- signif(daily_wy_aver2_589$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_589, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
589 SD
standard_dev_589 <- daily_wy_aver_589 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589 <- standard_dev_589 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_589 <- standard_dev_all_589 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_589 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.224944 |
| 1988 | 8.773281 |
| 1989 | 8.776693 |
| 1990 | 8.580289 |
| 1991 | 8.502502 |
| 1992 | 7.975611 |
| 1993 | 8.199081 |
| 1995 | 7.946318 |
| 1996 | 8.231469 |
| 1998 | 8.706109 |
| 1999 | 7.431979 |
| 2000 | 8.232988 |
| 2001 | 8.771085 |
| 2002 | 9.082003 |
| 2003 | 8.732769 |
| 2004 | 8.380419 |
| 2005 | 8.458900 |
| 2007 | 8.212611 |
| 2008 | 8.408220 |
| 2015 | 7.103666 |
| 2016 | 7.842413 |
| 2017 | 7.527960 |
| 2018 | 7.506884 |
| 2019 | 8.029199 |
| 2020 | 8.271099 |
| 2021 | 8.278625 |
| 2022 | 8.020852 |
ggplot(standard_dev_all_589, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 589 average temperatures for water years 2005-2021
MK & SS for 589 (non-corrected)
sd_mk_589 <- mk.test(standard_dev_all_589$sd_2)
print(sd_mk_589)##
## Mann-Kendall trend test
##
## data: standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -89.0000000 2301.0000000 -0.2535613sd_sens_589 <- sens.slope(standard_dev_all_589$sd_2)
print(sd_sens_589)##
## Sen's slope
##
## data: standard_dev_all_589$sd_2
## z = -1.8345, n = 27, p-value = 0.06658
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.043902404 0.001230005
## sample estimates:
## Sen's slope
## -0.02183399Corrected
589 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_589_ad <- snotel_589_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_589_ad <- yearly_wy_aver_589_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_589_ad <- daily_wy_aver_589_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_589_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_589_ad <-daily_wy_aver_589_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_589_ad$date_temp_ad <- signif(daily_wy_aver2_589_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_589_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
589 SS (corrected)
standard_dev_589_ad <- daily_wy_aver_589_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_589_ad <- standard_dev_589_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_589_ad <- standard_dev_all_589_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_589_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.671688 |
| 1988 | 8.200568 |
| 1989 | 8.222864 |
| 1990 | 7.999469 |
| 1991 | 7.962227 |
| 1992 | 7.420415 |
| 1993 | 7.639311 |
| 1995 | 7.378994 |
| 1996 | 7.641722 |
| 1998 | 8.086353 |
| 1999 | 6.914950 |
| 2000 | 7.633161 |
| 2001 | 8.156564 |
| 2002 | 8.440080 |
| 2003 | 8.103839 |
| 2004 | 7.827141 |
| 2005 | 7.781158 |
| 2007 | 8.212720 |
| 2008 | 8.408841 |
| 2015 | 7.104293 |
| 2016 | 7.842763 |
| 2017 | 7.527954 |
| 2018 | 7.507329 |
| 2019 | 8.029667 |
| 2020 | 8.271310 |
| 2021 | 8.278726 |
| 2022 | 8.021229 |
ggplot(standard_dev_all_589_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 589 average temperatures for water years 1986-2021
MK & SS 589 (corrected)
sd_mk_589_ad <- mk.test(standard_dev_all_589_ad$sd_2)
print(sd_mk_589_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 3.300000e+01 2.301000e+03 9.401709e-02sd_sens_589_ad <- sens.slope(standard_dev_all_589_ad$sd_2)
print(sd_sens_589_ad)##
## Sen's slope
##
## data: standard_dev_all_589_ad$sd_2
## z = 0.6671, n = 27, p-value = 0.5047
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01153608 0.02498422
## sample estimates:
## Sen's slope
## 0.008250776Middle Creek 624
Morrisey 6/26/2006
snotel_624 <- SNOTEL_san_juan_Area %>%
filter(site_id == "624")#str(snotel_624) # check the date, usually a character.
snotel_624$Date <- as.Date(snotel_624$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_624_clean <- snotel_624 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_624_clean <- snotel_624_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_624_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_624_clean <- snotel_624_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_624_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_624_cull_count <- snotel_624_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_624_cull_count## # A tibble: 43 x 2
## # Groups: waterYear [43]
## waterYear n
## <dbl> <int>
## 1 1980 59
## 2 1981 354
## 3 1982 1
## 4 1983 312
## 5 1984 360
## 6 1985 362
## 7 1986 364
## 8 1987 360
## 9 1988 340
## 10 1989 363
## # ... with 33 more rows# filtering for too few observations in a year
snotel_624_cull_count_days <- snotel_624_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_624_cull_count_days## # A tibble: 9 x 2
## # Groups: waterYear [9]
## waterYear n
## <dbl> <int>
## 1 1980 59
## 2 1982 1
## 3 1983 312
## 4 1988 340
## 5 1991 338
## 6 1993 330
## 7 1994 330
## 8 2021 338
## 9 2022 3481980, 1982, 1983, 1988, 1991, 1993, 1994, 2021, 2022 need to be culled.
snotel_624_clean_culled <- snotel_624_clean %>%
filter(waterYear != "1980" & waterYear != "1981" & waterYear != "1982" & waterYear != "1983" & waterYear != "1988" & waterYear != "1991" & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022")# & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_624_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_624_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)
dygraph(temp_624_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_624_clean_culled <- snotel_624_clean_culled %>%
filter(temperature_mean < 20)
temp_624_xts <- xts(snotel_624_clean_culled$temperature_mean, order.by = snotel_624_clean_culled$Date)
dygraph(temp_624_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Middle Creek 624 Morrisey 6/26/2006
snotel_624_adjusted <- snotel_624_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-06-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
624 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_624 <- snotel_624_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_624 <- yearly_wy_aver_624 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_624 <- daily_wy_aver_624 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_624$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_624 <-daily_wy_aver_624 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_624$date_temp <- signif(daily_wy_aver2_624$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_624, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
624 SD
standard_dev_624 <- daily_wy_aver_624 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624 <- standard_dev_624 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_624 <- standard_dev_all_624 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_624 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 8.239005 |
| 1985 | 8.207104 |
| 1986 | 7.191318 |
| 1987 | 7.871848 |
| 1989 | 8.454185 |
| 1990 | 8.200220 |
| 1992 | 6.438238 |
| 1995 | 7.565327 |
| 1996 | 7.835558 |
| 1997 | 7.996219 |
| 1998 | 8.415139 |
| 1999 | 7.174609 |
| 2000 | 7.839800 |
| 2001 | 8.355597 |
| 2002 | 8.492066 |
| 2003 | 8.150297 |
| 2004 | 7.912597 |
| 2005 | 7.915951 |
| 2006 | 7.936724 |
| 2007 | 7.640492 |
| 2008 | 8.021623 |
| 2009 | 7.171130 |
| 2010 | 8.259271 |
| 2011 | 7.932618 |
| 2012 | 7.791063 |
| 2013 | 8.259320 |
| 2014 | 7.372950 |
| 2015 | 6.806548 |
| 2016 | 7.548199 |
| 2017 | 7.108731 |
| 2018 | 7.036935 |
| 2019 | 7.675916 |
| 2020 | 7.741703 |
ggplot(standard_dev_all_624, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 624 average temperatures for water years 2005-2021
MK & SS for 624 (non-corrected)
sd_mk_624 <- mk.test(standard_dev_all_624$sd_2)
print(sd_mk_624)##
## Mann-Kendall trend test
##
## data: standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -128.0000000 4165.3333333 -0.2424242sd_sens_624 <- sens.slope(standard_dev_all_624$sd_2)
print(sd_sens_624)##
## Sen's slope
##
## data: standard_dev_all_624$sd_2
## z = -1.9678, n = 33, p-value = 0.04909
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.0333661824 -0.0003479011
## sample estimates:
## Sen's slope
## -0.01662145Corrected
624 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_624_ad <- snotel_624_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_624_ad <- yearly_wy_aver_624_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_624_ad <- daily_wy_aver_624_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_624_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_624_ad <-daily_wy_aver_624_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_624_ad$date_temp_ad <- signif(daily_wy_aver2_624_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_624_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
624 SS (corrected)
standard_dev_624_ad <- daily_wy_aver_624_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_624_ad <- standard_dev_624_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_624_ad <- standard_dev_all_624_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_624_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 7.738801 |
| 1985 | 7.713059 |
| 1986 | 6.725366 |
| 1987 | 7.374899 |
| 1989 | 7.942209 |
| 1990 | 7.687289 |
| 1992 | 6.080980 |
| 1995 | 7.080976 |
| 1996 | 7.336792 |
| 1997 | 7.519955 |
| 1998 | 7.885653 |
| 1999 | 6.718142 |
| 2000 | 7.329553 |
| 2001 | 7.833709 |
| 2002 | 7.944757 |
| 2003 | 7.617514 |
| 2004 | 7.442397 |
| 2005 | 7.390605 |
| 2006 | 7.304030 |
| 2007 | 7.639661 |
| 2008 | 8.021493 |
| 2009 | 7.170919 |
| 2010 | 8.259153 |
| 2011 | 7.932368 |
| 2012 | 7.790490 |
| 2013 | 8.259010 |
| 2014 | 7.372461 |
| 2015 | 6.806595 |
| 2016 | 7.547948 |
| 2017 | 7.107845 |
| 2018 | 7.036548 |
| 2019 | 7.675791 |
| 2020 | 7.741347 |
ggplot(standard_dev_all_624_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 624 average temperatures for water years 1986-2021
MK & SS 624 (corrected)
sd_mk_624_ad <- mk.test(standard_dev_all_624_ad$sd_2)
print(sd_mk_624_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 2.400000e+01 4.165333e+03 4.545455e-02sd_sens_624_ad <- sens.slope(standard_dev_all_624_ad$sd_2)
print(sd_sens_624_ad)##
## Sen's slope
##
## data: standard_dev_all_624_ad$sd_2
## z = 0.35637, n = 33, p-value = 0.7216
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01304163 0.02223675
## sample estimates:
## Sen's slope
## 0.003307866Molas Lake 632
NSCE-Morrisey, Bias-Original 10/2/2003
snotel_632 <- SNOTEL_san_juan_Area %>%
filter(site_id == "632")#str(snotel_632) # check the date, usually a character.
snotel_632$Date <- as.Date(snotel_632$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_632_clean <- snotel_632 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_632_clean <- snotel_632_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_632_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_632_clean <- snotel_632_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
#filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_632_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_632_cull_count <- snotel_632_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_632_cull_count## # A tibble: 37 x 2
## # Groups: waterYear [37]
## waterYear n
## <dbl> <int>
## 1 1986 54
## 2 1987 355
## 3 1988 362
## 4 1989 364
## 5 1990 365
## 6 1991 365
## 7 1992 366
## 8 1993 303
## 9 1994 294
## 10 1995 361
## # ... with 27 more rows# filtering for too few observations in a year
snotel_632_cull_count_days <- snotel_632_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_632_cull_count_days## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1986 54
## 2 1993 303
## 3 1994 295
## 4 2005 344
## 5 2006 304
## 6 2013 2481986, 1993, 1994, 2005, 2006, 2013 need to be culled.
snotel_632_clean_culled <- snotel_632_clean %>%
filter(waterYear != "1986" & waterYear != "1993" & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_632_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_632_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)
dygraph(temp_632_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_632_clean_culled <- snotel_632_clean_culled %>%
filter(temperature_mean < 20)
temp_632_xts <- xts(snotel_632_clean_culled$temperature_mean, order.by = snotel_632_clean_culled$Date)
dygraph(temp_632_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Molas Lake 632
NSCE-Morrisey, Bias-Original 10/2/2003
snotel_632_adjusted <- snotel_632_clean_culled %>%
mutate(temp_ad = if_else(Date < "2003-10-02", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
632 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_632 <- snotel_632_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_632 <- yearly_wy_aver_632 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_632 <- daily_wy_aver_632 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_632$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_632 <-daily_wy_aver_632 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_632$date_temp <- signif(daily_wy_aver2_632$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_632, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
632 SD
standard_dev_632 <- daily_wy_aver_632 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632 <- standard_dev_632 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_632 <- standard_dev_all_632 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_632 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.613234 |
| 1988 | 8.191410 |
| 1989 | 8.166070 |
| 1990 | 7.646926 |
| 1991 | 7.969691 |
| 1992 | 7.156808 |
| 1995 | 7.489888 |
| 1996 | 7.655841 |
| 1997 | 7.698012 |
| 1998 | 8.157837 |
| 1999 | 6.885991 |
| 2000 | 7.657554 |
| 2001 | 6.788385 |
| 2002 | 8.254532 |
| 2003 | 7.907020 |
| 2004 | 7.308264 |
| 2007 | 7.607147 |
| 2008 | 7.964946 |
| 2009 | 6.957483 |
| 2010 | 7.983696 |
| 2011 | 7.701623 |
| 2012 | 7.490054 |
| 2014 | 7.236529 |
| 2015 | 6.581104 |
| 2016 | 7.409875 |
| 2017 | 7.001441 |
| 2018 | 6.814621 |
| 2019 | 7.582515 |
| 2020 | 7.766470 |
| 2021 | 7.792571 |
| 2022 | 7.357711 |
ggplot(standard_dev_all_632, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 632 average temperatures for water years 2005-2021
MK & SS for 632 (non-corrected)
sd_mk_632 <- mk.test(standard_dev_all_632$sd_2)
print(sd_mk_632)##
## Mann-Kendall trend test
##
## data: standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -97.0000000 3461.6666667 -0.2086022sd_sens_632 <- sens.slope(standard_dev_all_632$sd_2)
print(sd_sens_632)##
## Sen's slope
##
## data: standard_dev_all_632$sd_2
## z = -1.6317, n = 31, p-value = 0.1028
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.034372524 0.003422909
## sample estimates:
## Sen's slope
## -0.01583493Corrected
632 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_632_ad <- snotel_632_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_632_ad <- yearly_wy_aver_632_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_632_ad <- daily_wy_aver_632_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_632_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_632_ad <-daily_wy_aver_632_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_632_ad$date_temp_ad <- signif(daily_wy_aver2_632_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_632_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
632 SS (corrected)
standard_dev_632_ad <- daily_wy_aver_632_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_632_ad <- standard_dev_632_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_632_ad <- standard_dev_all_632_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_632_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.106416 |
| 1988 | 7.678134 |
| 1989 | 7.662133 |
| 1990 | 7.141729 |
| 1991 | 7.492618 |
| 1992 | 6.671512 |
| 1995 | 6.957618 |
| 1996 | 7.110780 |
| 1997 | 7.178235 |
| 1998 | 7.586041 |
| 1999 | 6.397699 |
| 2000 | 7.103150 |
| 2001 | 6.358660 |
| 2002 | 7.672831 |
| 2003 | 7.338249 |
| 2004 | 7.309361 |
| 2007 | 7.606845 |
| 2008 | 7.965266 |
| 2009 | 6.957636 |
| 2010 | 7.984148 |
| 2011 | 7.701741 |
| 2012 | 7.490245 |
| 2014 | 7.236308 |
| 2015 | 6.581364 |
| 2016 | 7.409996 |
| 2017 | 7.000722 |
| 2018 | 6.815246 |
| 2019 | 7.583475 |
| 2020 | 7.766635 |
| 2021 | 7.793272 |
| 2022 | 7.358050 |
ggplot(standard_dev_all_632_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 632 average temperatures for water years 1986-2021
MK & SS 632 (corrected)
sd_mk_632_ad <- mk.test(standard_dev_all_632_ad$sd_2)
print(sd_mk_632_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 57.0000000 3461.6666667 0.1225806sd_sens_632_ad <- sens.slope(standard_dev_all_632_ad$sd_2)
print(sd_sens_632_ad)##
## Sen's slope
##
## data: standard_dev_all_632_ad$sd_2
## z = 0.9518, n = 31, p-value = 0.3412
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01086011 0.02505935
## sample estimates:
## Sen's slope
## 0.008387789Red Mountain Pass 713
Morrisey 8/18/2004
snotel_713 <- SNOTEL_san_juan_Area %>%
filter(site_id == "713")#str(snotel_713) # check the date, usually a character.
snotel_713$Date <- as.Date(snotel_713$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_713_clean <- snotel_713 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_713_clean <- snotel_713_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_713_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_713_clean <- snotel_713_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_713_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_713_cull_count <- snotel_713_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_713_cull_count## # A tibble: 42 x 2
## # Groups: waterYear [42]
## waterYear n
## <dbl> <int>
## 1 1981 189
## 2 1982 287
## 3 1983 319
## 4 1984 360
## 5 1985 362
## 6 1986 364
## 7 1987 263
## 8 1988 282
## 9 1989 228
## 10 1990 365
## # ... with 32 more rows# filtering for too few observations in a year
snotel_713_cull_count_days <- snotel_713_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_713_cull_count_days## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1981 189
## 2 1982 287
## 3 1983 319
## 4 1987 315
## 5 1989 263
## 6 1994 3081981, 1982, 1983, 1987, 1989, 1994 need to be culled. 1984 & 1985 are too low.
snotel_713_clean_culled <- snotel_713_clean %>%
filter(waterYear > "1985" & waterYear != "1987" & waterYear != "1989" & waterYear != "1994")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_713_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_713_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_713_xts <- xts(snotel_713_clean_culled$temperature_mean, order.by = snotel_713_clean_culled$Date)
dygraph(temp_713_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_713_clean_culled <- snotel_713_clean_culled %>%
filter(temperature_mean < 19.3)
temp_713_xts <- xts(snotel_713_clean_culled$temperature_mean, order.by = snotel_713_clean_culled$Date)
dygraph(temp_713_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Red Mountain Pass 713
Morrisey 8/18/2004
snotel_713_adjusted <- snotel_713_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-08-18", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
713 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_713 <- snotel_713_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_713 <- yearly_wy_aver_713 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_713 <- daily_wy_aver_713 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_713$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_713 <-daily_wy_aver_713 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_713$date_temp <- signif(daily_wy_aver2_713$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_713, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
713 SD
standard_dev_713 <- daily_wy_aver_713 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_713 <- standard_dev_713 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_713 <- standard_dev_all_713 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_713 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1986 | 7.022683 |
| 1988 | 9.344095 |
| 1990 | 7.761743 |
| 1991 | 7.768555 |
| 1992 | 7.212083 |
| 1993 | 7.716280 |
| 1995 | 7.336781 |
| 1996 | 7.731429 |
| 1997 | 7.947757 |
| 1998 | 8.198019 |
| 1999 | 6.921528 |
| 2000 | 7.680921 |
| 2001 | 8.089853 |
| 2002 | 8.395619 |
| 2003 | 7.952564 |
| 2004 | 7.777488 |
| 2005 | 7.151965 |
| 2006 | 7.312223 |
| 2007 | 7.636775 |
| 2008 | 7.908986 |
| 2009 | 6.940732 |
| 2010 | 7.877650 |
| 2011 | 7.730288 |
| 2012 | 7.531887 |
| 2013 | 8.195738 |
| 2014 | 7.284701 |
| 2015 | 6.638343 |
| 2016 | 7.380700 |
| 2017 | 7.097935 |
| 2018 | 7.002862 |
| 2019 | 7.626598 |
| 2020 | 7.810018 |
| 2021 | 7.817233 |
| 2022 | 7.339213 |
ggplot(standard_dev_all_713, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 713 average temperatures for water years 2005-2021
MK & SS for 713 (non-corrected)
sd_mk_713 <- mk.test(standard_dev_all_713$sd_2)
print(sd_mk_713)##
## Mann-Kendall trend test
##
## data: standard_dev_all_713$sd_2
## z = -1.1563, n = 34, p-value = 0.2476
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -79.00000 4550.33333 -0.14082sd_sens_713 <- sens.slope(standard_dev_all_713$sd_2)
print(sd_sens_713)##
## Sen's slope
##
## data: standard_dev_all_713$sd_2
## z = -1.1563, n = 34, p-value = 0.2476
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.028289222 0.006454843
## sample estimates:
## Sen's slope
## -0.008980292Corrected
713 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_713_ad <- snotel_713_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_713_ad <- yearly_wy_aver_713_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_713_ad <- daily_wy_aver_713_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_713_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_713_ad <-daily_wy_aver_713_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_713_ad$date_temp_ad <- signif(daily_wy_aver2_713_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_713_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
713 SS (corrected)
standard_dev_713_ad <- daily_wy_aver_713_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_713_ad <- standard_dev_713_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_713_ad <- standard_dev_all_713_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_713_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1986 | 6.564707 |
| 1988 | 8.872885 |
| 1990 | 7.267850 |
| 1991 | 7.311859 |
| 1992 | 6.739882 |
| 1993 | 7.221869 |
| 1995 | 6.846024 |
| 1996 | 7.223659 |
| 1997 | 7.456352 |
| 1998 | 7.661532 |
| 1999 | 6.466454 |
| 2000 | 7.159552 |
| 2001 | 7.557097 |
| 2002 | 7.840372 |
| 2003 | 7.411227 |
| 2004 | 7.238021 |
| 2005 | 7.151835 |
| 2006 | 7.312105 |
| 2007 | 7.636625 |
| 2008 | 7.909548 |
| 2009 | 6.940894 |
| 2010 | 7.878182 |
| 2011 | 7.730889 |
| 2012 | 7.532385 |
| 2013 | 8.196125 |
| 2014 | 7.284711 |
| 2015 | 6.639012 |
| 2016 | 7.380681 |
| 2017 | 7.097919 |
| 2018 | 7.003100 |
| 2019 | 7.627224 |
| 2020 | 7.810395 |
| 2021 | 7.817455 |
| 2022 | 7.339782 |
ggplot(standard_dev_all_713_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 713 average temperatures for water years 1986-2021
MK & SS 713 (corrected)
sd_mk_713_ad <- mk.test(standard_dev_all_713_ad$sd_2)
print(sd_mk_713_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_713_ad$sd_2
## z = 1.4231, n = 34, p-value = 0.1547
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 97.0000000 4550.3333333 0.1729055sd_sens_713_ad <- sens.slope(standard_dev_all_713_ad$sd_2)
print(sd_sens_713_ad)##
## Sen's slope
##
## data: standard_dev_all_713_ad$sd_2
## z = 1.4231, n = 34, p-value = 0.1547
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.005520634 0.030221234
## sample estimates:
## Sen's slope
## 0.01194811Scotch Creek 739
Morrisey 6/15/2004
snotel_739 <- SNOTEL_san_juan_Area %>%
filter(site_id == "739")#str(snotel_739) # check the date, usually a character.
snotel_739$Date <- as.Date(snotel_739$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_739_clean <- snotel_739 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_739_clean <- snotel_739_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_739_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_739_clean <- snotel_739_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_739_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
snotel_739_cull_count <- snotel_739_clean %>%
filter(temperature_min > -40) %>%
count(waterYear)
snotel_739_cull_count## # A tibble: 36 x 2
## # Groups: waterYear [36]
## waterYear n
## <dbl> <int>
## 1 1987 345
## 2 1988 273
## 3 1989 363
## 4 1990 365
## 5 1991 365
## 6 1992 365
## 7 1993 342
## 8 1994 345
## 9 1995 362
## 10 1996 365
## # ... with 26 more rows# filtering for too few observations in a year
snotel_739_cull_count_days <- snotel_739_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_739_cull_count_days## # A tibble: 6 x 2
## # Groups: waterYear [6]
## waterYear n
## <dbl> <int>
## 1 1987 345
## 2 1988 273
## 3 1993 343
## 4 1994 345
## 5 2016 320
## 6 2017 2491987, 1988, 1993, 1994, 2016, 2017 need to be culled.
snotel_739_clean_culled <- snotel_739_clean %>%
filter(waterYear != "1987" & waterYear != "1988" & waterYear != "1993" & waterYear != "1994" & waterYear != "2016" & waterYear != "2017")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_739_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_739_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_739_xts <- xts(snotel_739_clean_culled$temperature_mean, order.by = snotel_739_clean_culled$Date)
dygraph(temp_739_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") snotel_739_clean_culled <- snotel_739_clean_culled %>%
filter(temperature_mean < 19.3)
temp_739_xts <- xts(snotel_739_clean_culled$temperature_mean, order.by = snotel_739_clean_culled$Date)
dygraph(temp_739_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Scotch Creek 739
Morrisey 6/15/2004
snotel_739_adjusted <- snotel_739_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-06-15", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
739 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_739 <- snotel_739_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_739 <- yearly_wy_aver_739 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_739 <- daily_wy_aver_739 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_739$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_739 <-daily_wy_aver_739 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_739$date_temp <- signif(daily_wy_aver2_739$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_739, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
739 SD
standard_dev_739 <- daily_wy_aver_739 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_739 <- standard_dev_739 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_739 <- standard_dev_all_739 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_739 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1989 | 7.327097 |
| 1990 | 6.021389 |
| 1991 | 5.989983 |
| 1992 | 4.863013 |
| 1995 | 7.600315 |
| 1996 | 7.840086 |
| 1997 | 7.952598 |
| 1998 | 8.393336 |
| 1999 | 7.124567 |
| 2000 | 7.697346 |
| 2001 | 8.108824 |
| 2002 | 8.727256 |
| 2003 | 7.977395 |
| 2004 | 7.894210 |
| 2005 | 7.144352 |
| 2006 | 7.250093 |
| 2007 | 7.615681 |
| 2008 | 7.916198 |
| 2009 | 6.788880 |
| 2010 | 7.880288 |
| 2011 | 7.703849 |
| 2012 | 7.458561 |
| 2013 | 8.171355 |
| 2014 | 7.223604 |
| 2015 | 6.545241 |
| 2018 | 6.908918 |
| 2019 | 7.558347 |
| 2020 | 7.673725 |
| 2021 | 7.658883 |
| 2022 | 7.427996 |
ggplot(standard_dev_all_739, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 739 average temperatures for water years 2005-2021
MK & SS for 739 (non-corrected)
sd_mk_739 <- mk.test(standard_dev_all_739$sd_2)
print(sd_mk_739)##
## Mann-Kendall trend test
##
## data: standard_dev_all_739$sd_2
## z = 0, n = 30, p-value = 1
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 1.000000e+00 3.141667e+03 2.298851e-03sd_sens_739 <- sens.slope(standard_dev_all_739$sd_2)
print(sd_sens_739)##
## Sen's slope
##
## data: standard_dev_all_739$sd_2
## z = 0, n = 30, p-value = 1
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.02559406 0.03474477
## sample estimates:
## Sen's slope
## 0.0005911373Corrected
739 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_739_ad <- snotel_739_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_739_ad <- yearly_wy_aver_739_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_739_ad <- daily_wy_aver_739_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_739_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_739_ad <-daily_wy_aver_739_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_739_ad$date_temp_ad <- signif(daily_wy_aver2_739_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_739_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
739 SS (corrected)
standard_dev_739_ad <- daily_wy_aver_739_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_739_ad <- standard_dev_739_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_739_ad <- standard_dev_all_739_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_739_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1989 | 6.882254 |
| 1990 | 5.641380 |
| 1991 | 5.682802 |
| 1992 | 4.573743 |
| 1995 | 7.029026 |
| 1996 | 7.246834 |
| 1997 | 7.374205 |
| 1998 | 7.773434 |
| 1999 | 6.583515 |
| 2000 | 7.110733 |
| 2001 | 7.493196 |
| 2002 | 8.076937 |
| 2003 | 7.362576 |
| 2004 | 7.247292 |
| 2005 | 7.144310 |
| 2006 | 7.250207 |
| 2007 | 7.615515 |
| 2008 | 7.916383 |
| 2009 | 6.789140 |
| 2010 | 7.880400 |
| 2011 | 7.703781 |
| 2012 | 7.458914 |
| 2013 | 8.171585 |
| 2014 | 7.223615 |
| 2015 | 6.545250 |
| 2018 | 6.908860 |
| 2019 | 7.558520 |
| 2020 | 7.674026 |
| 2021 | 7.658936 |
| 2022 | 7.428362 |
ggplot(standard_dev_all_739_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 739 average temperatures for water years 1986-2021
MK & SS 739 (corrected)
sd_mk_739_ad <- mk.test(standard_dev_all_739_ad$sd_2)
print(sd_mk_739_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_739_ad$sd_2
## z = 2.3907, n = 30, p-value = 0.01682
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 135.0000000 3141.6666667 0.3103448sd_sens_739_ad <- sens.slope(standard_dev_all_739_ad$sd_2)
print(sd_sens_739_ad)##
## Sen's slope
##
## data: standard_dev_all_739_ad$sd_2
## z = 2.3907, n = 30, p-value = 0.01682
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## 0.00806300 0.05860595
## sample estimates:
## Sen's slope
## 0.02802846Slumgullion 762
Morrisey 6/26/2006
snotel_762 <- SNOTEL_san_juan_Area %>%
filter(site_id == "762")#str(snotel_762) # check the date, usually a character.
snotel_762$Date <- as.Date(snotel_762$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_762_clean <- snotel_762 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_762_clean <- snotel_762_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_762_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_762_clean <- snotel_762_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_762_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_762_cull_count <- snotel_762_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_762_cull_count
# filtering for too few observations in a year
snotel_762_cull_count_days <- snotel_762_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_762_cull_count_days## # A tibble: 5 x 2
## # Groups: waterYear [5]
## waterYear n
## <dbl> <int>
## 1 1980 43
## 2 1982 1
## 3 1983 319
## 4 1994 331
## 5 2003 3371980, 1982, 1994, 2003 need to be culled. 1995 maximums are less than the minimums.
snotel_762_clean_culled <- snotel_762_clean %>%
filter(waterYear != "1980" & waterYear != "1981" & waterYear != "1982" & waterYear != "1983" & waterYear != "1994" & waterYear != "1995" & waterYear != "2003")# & waterYear != "2017" & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_762_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_762_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_762_xts <- xts(snotel_762_clean_culled$temperature_mean, order.by = snotel_762_clean_culled$Date)
dygraph(temp_762_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_762_clean_culled <- snotel_762_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_762_xts <- xts(snotel_762_clean_culled$temperature_mean, order.by = snotel_762_clean_culled$Date)
dygraph(temp_762_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Slumgullion 762
Morrisey 6/26/2006
snotel_762_adjusted <- snotel_762_clean_culled %>%
mutate(temp_ad = if_else(Date < "2006-06-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
762 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_762 <- snotel_762_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_762 <- yearly_wy_aver_762 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_762 <- daily_wy_aver_762 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_762$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_762 <-daily_wy_aver_762 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_762$date_temp <- signif(daily_wy_aver2_762$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_762, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
762 SD
standard_dev_762 <- daily_wy_aver_762 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_762 <- standard_dev_762 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_762 <- standard_dev_all_762 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_762 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 8.474614 |
| 1985 | 8.395416 |
| 1986 | 7.224771 |
| 1987 | 8.036169 |
| 1988 | 8.470217 |
| 1989 | 8.556810 |
| 1990 | 8.104917 |
| 1991 | 8.129177 |
| 1992 | 7.378858 |
| 1993 | 7.881362 |
| 1996 | 8.151725 |
| 1997 | 8.002402 |
| 1998 | 8.492764 |
| 1999 | 7.140044 |
| 2000 | 7.859908 |
| 2001 | 8.456394 |
| 2002 | 8.777834 |
| 2004 | 8.172096 |
| 2005 | 8.157803 |
| 2006 | 8.479381 |
| 2007 | 7.832136 |
| 2008 | 8.378519 |
| 2009 | 7.328609 |
| 2010 | 8.218219 |
| 2011 | 8.186856 |
| 2012 | 7.846460 |
| 2013 | 8.493891 |
| 2014 | 7.772861 |
| 2015 | 7.080609 |
| 2016 | 7.803357 |
| 2017 | 7.532567 |
| 2018 | 7.423620 |
| 2019 | 8.010211 |
| 2020 | 8.142119 |
| 2021 | 8.245368 |
| 2022 | 7.614681 |
ggplot(standard_dev_all_762, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 762 average temperatures for water years 2005-2021
MK & SS for 762 (non-corrected)
sd_mk_762 <- mk.test(standard_dev_all_762$sd_2)
print(sd_mk_762)##
## Mann-Kendall trend test
##
## data: standard_dev_all_762$sd_2
## z = -1.4574, n = 36, p-value = 0.145
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -108.0000000 5390.0000000 -0.1714286sd_sens_762 <- sens.slope(standard_dev_all_762$sd_2)
print(sd_sens_762)##
## Sen's slope
##
## data: standard_dev_all_762$sd_2
## z = -1.4574, n = 36, p-value = 0.145
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.024958596 0.002910121
## sample estimates:
## Sen's slope
## -0.01112266Corrected
762 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_762_ad <- snotel_762_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_762_ad <- yearly_wy_aver_762_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_762_ad <- daily_wy_aver_762_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_762_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_762_ad <-daily_wy_aver_762_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_762_ad$date_temp_ad <- signif(daily_wy_aver2_762_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_762_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
762 SS (corrected)
standard_dev_762_ad <- daily_wy_aver_762_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_762_ad <- standard_dev_762_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_762_ad <- standard_dev_all_762_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_762_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1984 | 7.970791 |
| 1985 | 7.888311 |
| 1986 | 6.764665 |
| 1987 | 7.528970 |
| 1988 | 7.965877 |
| 1989 | 8.053090 |
| 1990 | 7.604786 |
| 1991 | 7.665906 |
| 1992 | 6.904408 |
| 1993 | 7.400370 |
| 1996 | 7.629529 |
| 1997 | 7.511155 |
| 1998 | 7.956453 |
| 1999 | 6.688385 |
| 2000 | 7.353371 |
| 2001 | 7.919104 |
| 2002 | 8.206929 |
| 2004 | 7.685091 |
| 2005 | 7.619040 |
| 2006 | 7.823927 |
| 2007 | 7.831876 |
| 2008 | 8.379138 |
| 2009 | 7.329078 |
| 2010 | 8.218734 |
| 2011 | 8.187217 |
| 2012 | 7.846411 |
| 2013 | 8.493922 |
| 2014 | 7.773398 |
| 2015 | 7.081459 |
| 2016 | 7.804132 |
| 2017 | 7.532344 |
| 2018 | 7.424275 |
| 2019 | 8.010639 |
| 2020 | 8.142599 |
| 2021 | 8.245690 |
| 2022 | 7.615524 |
ggplot(standard_dev_all_762_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 762 average temperatures for water years 1986-2021
MK & SS 762 (corrected)
sd_mk_762_ad <- mk.test(standard_dev_all_762_ad$sd_2)
print(sd_mk_762_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_762_ad$sd_2
## z = 1.2123, n = 36, p-value = 0.2254
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 90.0000000 5390.0000000 0.1428571sd_sens_762_ad <- sens.slope(standard_dev_all_762_ad$sd_2)
print(sd_sens_762_ad)##
## Sen's slope
##
## data: standard_dev_all_762_ad$sd_2
## z = 1.2123, n = 36, p-value = 0.2254
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.00495755 0.02329061
## sample estimates:
## Sen's slope
## 0.009103679Spud Mountain 780
Morrisey 6/28/2004
snotel_780 <- SNOTEL_san_juan_Area %>%
filter(site_id == "780")#str(snotel_780) # check the date, usually a character.
snotel_780$Date <- as.Date(snotel_780$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_780_clean <- snotel_780 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_780_clean <- snotel_780_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_780_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_780_clean <- snotel_780_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_780_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_780_cull_count <- snotel_780_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_780_cull_count
# filtering for too few observations in a year
snotel_780_cull_count_days <- snotel_780_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_780_cull_count_days## # A tibble: 8 x 2
## # Groups: waterYear [8]
## waterYear n
## <dbl> <int>
## 1 1987 315
## 2 1994 336
## 3 1998 309
## 4 2009 342
## 5 2011 337
## 6 2013 336
## 7 2016 331
## 8 2017 3241987, 1994, 1998, 2009, 2011, 2013, 2016, 2017 need to be culled.
snotel_780_clean_culled <- snotel_780_clean %>%
filter(waterYear != "1987" & waterYear != "1994" & waterYear != "1998" & waterYear != "2009" & waterYear != "2011" & waterYear != "2013" & waterYear != "2016" & waterYear != "2017")# & waterYear != "1993" & waterYear != "1994" & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_780_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_780_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_780_xts <- xts(snotel_780_clean_culled$temperature_mean, order.by = snotel_780_clean_culled$Date)
dygraph(temp_780_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_780_clean_culled <- snotel_780_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_780_xts <- xts(snotel_780_clean_culled$temperature_mean, order.by = snotel_780_clean_culled$Date)
dygraph(temp_780_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Spud Mountain 780
Morrisey 6/28/2004
snotel_780_adjusted <- snotel_780_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-06-28", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
780 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_780 <- snotel_780_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_780 <- yearly_wy_aver_780 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_780 <- daily_wy_aver_780 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_780$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_780 <-daily_wy_aver_780 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_780$date_temp <- signif(daily_wy_aver2_780$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_780, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
780 SD
standard_dev_780 <- daily_wy_aver_780 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_780 <- standard_dev_780 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_780 <- standard_dev_all_780 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_780 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 8.083874 |
| 1989 | 8.173292 |
| 1990 | 7.631378 |
| 1991 | 7.796215 |
| 1992 | 7.011676 |
| 1993 | 7.893616 |
| 1995 | 7.478626 |
| 1996 | 7.773665 |
| 1997 | 7.856295 |
| 1999 | 6.828980 |
| 2000 | 7.722333 |
| 2001 | 8.147987 |
| 2002 | 8.545717 |
| 2003 | 7.955225 |
| 2004 | 7.944675 |
| 2005 | 7.146399 |
| 2006 | 7.088438 |
| 2007 | 7.424508 |
| 2008 | 7.859459 |
| 2010 | 7.733705 |
| 2012 | 7.461301 |
| 2014 | 7.143839 |
| 2015 | 6.533563 |
| 2018 | 6.864656 |
| 2019 | 7.580081 |
| 2020 | 7.727076 |
| 2021 | 7.724978 |
| 2022 | 7.146411 |
ggplot(standard_dev_all_780, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 780 average temperatures for water years 2005-2021
MK & SS for 780 (non-corrected)
sd_mk_780 <- mk.test(standard_dev_all_780$sd_2)
print(sd_mk_780)##
## Mann-Kendall trend test
##
## data: standard_dev_all_780$sd_2
## z = -1.9954, n = 28, p-value = 0.046
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -102.0000000 2562.0000000 -0.2698413sd_sens_780 <- sens.slope(standard_dev_all_780$sd_2)
print(sd_sens_780)##
## Sen's slope
##
## data: standard_dev_all_780$sd_2
## z = -1.9954, n = 28, p-value = 0.046
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.0394954014 -0.0004267439
## sample estimates:
## Sen's slope
## -0.01908281Corrected
780 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_780_ad <- snotel_780_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_780_ad <- yearly_wy_aver_780_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_780_ad <- daily_wy_aver_780_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_780_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_780_ad <-daily_wy_aver_780_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_780_ad$date_temp_ad <- signif(daily_wy_aver2_780_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_780_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
780 SS (corrected)
standard_dev_780_ad <- daily_wy_aver_780_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_780_ad <- standard_dev_780_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_780_ad <- standard_dev_all_780_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_780_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 7.530892 |
| 1989 | 7.629513 |
| 1990 | 7.085085 |
| 1991 | 7.275264 |
| 1992 | 6.494536 |
| 1993 | 7.336765 |
| 1995 | 6.928211 |
| 1996 | 7.209777 |
| 1997 | 7.307706 |
| 1999 | 6.331029 |
| 2000 | 7.141268 |
| 2001 | 7.555861 |
| 2002 | 7.925084 |
| 2003 | 7.356618 |
| 2004 | 7.299749 |
| 2005 | 7.146735 |
| 2006 | 7.089360 |
| 2007 | 7.425654 |
| 2008 | 7.860535 |
| 2010 | 7.734756 |
| 2012 | 7.462985 |
| 2014 | 7.144698 |
| 2015 | 6.534943 |
| 2018 | 6.866003 |
| 2019 | 7.581136 |
| 2020 | 7.728271 |
| 2021 | 7.725902 |
| 2022 | 7.147773 |
ggplot(standard_dev_all_780_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 780 average temperatures for water years 1986-2021
MK & SS 780 (corrected)
sd_mk_780_ad <- mk.test(standard_dev_all_780_ad$sd_2)
print(sd_mk_780_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_780_ad$sd_2
## z = 0.84953, n = 28, p-value = 0.3956
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 44.0000000 2562.0000000 0.1164021sd_sens_780_ad <- sens.slope(standard_dev_all_780_ad$sd_2)
print(sd_sens_780_ad)##
## Sen's slope
##
## data: standard_dev_all_780_ad$sd_2
## z = 0.84953, n = 28, p-value = 0.3956
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.01274116 0.02720599
## sample estimates:
## Sen's slope
## 0.007453894Stump Lakes 797
Morrisey 7/22/2005
snotel_797 <- SNOTEL_san_juan_Area %>%
filter(site_id == "797")#str(snotel_797) # check the date, usually a character.
snotel_797$Date <- as.Date(snotel_797$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_797_clean <- snotel_797 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_797_clean <- snotel_797_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_797_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_797_clean <- snotel_797_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 40)ggplot(snotel_797_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_797_cull_count <- snotel_797_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_797_cull_count
# filtering for too few observations in a year
snotel_797_cull_count_days <- snotel_797_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_797_cull_count_days## # A tibble: 5 x 2
## # Groups: waterYear [5]
## waterYear n
## <dbl> <int>
## 1 1987 267
## 2 1993 320
## 3 1994 274
## 4 2005 329
## 5 2006 3271987, 1993, 1994, 2005, 2006 need to be culled.
snotel_797_clean_culled <- snotel_797_clean %>%
filter(waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "1993" & waterYear != "1987")# & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_797_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_797_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_797_xts <- xts(snotel_797_clean_culled$temperature_mean, order.by = snotel_797_clean_culled$Date)
dygraph(temp_797_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_797_clean_culled <- snotel_797_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_797_xts <- xts(snotel_797_clean_culled$temperature_mean, order.by = snotel_797_clean_culled$Date)
dygraph(temp_797_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Stump Lakes 797
Morrisey 7/22/2005
snotel_797_adjusted <- snotel_797_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-07-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
797 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_797 <- snotel_797_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_797 <- yearly_wy_aver_797 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_797 <- daily_wy_aver_797 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_797$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_797 <-daily_wy_aver_797 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_797$date_temp <- signif(daily_wy_aver2_797$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_797, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
797 SD
standard_dev_797 <- daily_wy_aver_797 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_797 <- standard_dev_797 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_797 <- standard_dev_all_797 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_797 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 8.133643 |
| 1989 | 10.300128 |
| 1990 | 7.808751 |
| 1991 | 9.471092 |
| 1992 | 7.312520 |
| 1995 | 7.694636 |
| 1996 | 7.953496 |
| 1997 | 8.069085 |
| 1998 | 8.367447 |
| 1999 | 7.015052 |
| 2000 | 7.935172 |
| 2001 | 8.545423 |
| 2002 | 9.270169 |
| 2003 | 8.084230 |
| 2004 | 9.099232 |
| 2007 | 7.603452 |
| 2008 | 8.008843 |
| 2009 | 7.184718 |
| 2010 | 7.913867 |
| 2011 | 7.886632 |
| 2012 | 7.696409 |
| 2013 | 8.096842 |
| 2014 | 7.327713 |
| 2015 | 6.737528 |
| 2016 | 7.565177 |
| 2017 | 7.316786 |
| 2018 | 7.076498 |
| 2019 | 7.800277 |
| 2020 | 7.872608 |
| 2021 | 7.863889 |
| 2022 | 7.341735 |
ggplot(standard_dev_all_797, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 797 average temperatures for water years 2005-2021
MK & SS for 797 (non-corrected)
sd_mk_797 <- mk.test(standard_dev_all_797$sd_2)
print(sd_mk_797)##
## Mann-Kendall trend test
##
## data: standard_dev_all_797$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -151.0000000 3461.6666667 -0.3247312sd_sens_797 <- sens.slope(standard_dev_all_797$sd_2)
print(sd_sens_797)##
## Sen's slope
##
## data: standard_dev_all_797$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.055919889 -0.005393295
## sample estimates:
## Sen's slope
## -0.02911914Corrected
797 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_797_ad <- snotel_797_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_797_ad <- yearly_wy_aver_797_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_797_ad <- daily_wy_aver_797_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_797_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_797_ad <-daily_wy_aver_797_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_797_ad$date_temp_ad <- signif(daily_wy_aver2_797_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_797_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
797 SS (corrected)
standard_dev_797_ad <- daily_wy_aver_797_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_797_ad <- standard_dev_797_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_797_ad <- standard_dev_all_797_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_797_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1988 | 7.629765 |
| 1989 | 9.567031 |
| 1990 | 7.287673 |
| 1991 | 8.816683 |
| 1992 | 6.729102 |
| 1995 | 7.165683 |
| 1996 | 7.422437 |
| 1997 | 7.545973 |
| 1998 | 7.813938 |
| 1999 | 6.542177 |
| 2000 | 7.377392 |
| 2001 | 7.978790 |
| 2002 | 8.614161 |
| 2003 | 7.487087 |
| 2004 | 8.402698 |
| 2007 | 7.603448 |
| 2008 | 8.008551 |
| 2009 | 7.184680 |
| 2010 | 7.914300 |
| 2011 | 7.886697 |
| 2012 | 7.696830 |
| 2013 | 8.096444 |
| 2014 | 7.327715 |
| 2015 | 6.737458 |
| 2016 | 7.565155 |
| 2017 | 7.316299 |
| 2018 | 7.076457 |
| 2019 | 7.800532 |
| 2020 | 7.873011 |
| 2021 | 7.863666 |
| 2022 | 7.341402 |
ggplot(standard_dev_all_797_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 797 average temperatures for water years 1986-2021
MK & SS 797 (corrected)
sd_mk_797_ad <- mk.test(standard_dev_all_797_ad$sd_2)
print(sd_mk_797_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_797_ad$sd_2
## z = -0.30594, n = 31, p-value = 0.7597
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -19.00000000 3461.66666667 -0.04086022sd_sens_797_ad <- sens.slope(standard_dev_all_797_ad$sd_2)
print(sd_sens_797_ad)##
## Sen's slope
##
## data: standard_dev_all_797_ad$sd_2
## z = -0.30594, n = 31, p-value = 0.7597
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.03291859 0.02011237
## sample estimates:
## Sen's slope
## -0.00412891Upper Rio Grande 839
Oyler -> Morrisey 5/26/2004
snotel_839 <- SNOTEL_san_juan_Area %>%
filter(site_id == "839")#str(snotel_839) # check the date, usually a character.
snotel_839$Date <- as.Date(snotel_839$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_839_clean <- snotel_839 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_839_clean <- snotel_839_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_839_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_839_clean <- snotel_839_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_839_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_839_cull_count <- snotel_839_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_839_cull_count
# filtering for too few observations in a year
snotel_839_cull_count_days <- snotel_839_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_839_cull_count_days## # A tibble: 5 x 2
## # Groups: waterYear [5]
## waterYear n
## <dbl> <int>
## 1 1994 338
## 2 2004 330
## 3 2012 339
## 4 2014 313
## 5 2016 3471994, 2004, 2012, 2014, 2016 need to be culled.
snotel_839_clean_culled <- snotel_839_clean %>%
filter(waterYear != "1994" & waterYear != "2004" & waterYear != "2012" & waterYear != "2014" & waterYear != "2016")# & waterYear != "2021" & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_839_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_839_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_839_xts <- xts(snotel_839_clean_culled$temperature_mean, order.by = snotel_839_clean_culled$Date)
dygraph(temp_839_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_839_clean_culled <- snotel_839_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_839_xts <- xts(snotel_839_clean_culled$temperature_mean, order.by = snotel_839_clean_culled$Date)
dygraph(temp_839_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Upper Rio Grande 839
Oyler -> Morrisey 5/26/2004
snotel_839_adjusted <- snotel_839_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-05-26", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
839 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_839 <- snotel_839_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_839 <- yearly_wy_aver_839 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_839 <- daily_wy_aver_839 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_839$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_839 <-daily_wy_aver_839 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_839$date_temp <- signif(daily_wy_aver2_839$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_839, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
839 SD
standard_dev_839 <- daily_wy_aver_839 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_839 <- standard_dev_839 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_839 <- standard_dev_all_839 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_839 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 9.000013 |
| 1988 | 9.231875 |
| 1989 | 9.553651 |
| 1990 | 8.539684 |
| 1991 | 9.095321 |
| 1992 | 8.840589 |
| 1993 | 9.034837 |
| 1995 | 8.477482 |
| 1996 | 8.374844 |
| 1997 | 8.635328 |
| 1998 | 9.151866 |
| 1999 | 7.826807 |
| 2000 | 8.411485 |
| 2001 | 9.307564 |
| 2002 | 9.394376 |
| 2003 | 8.653297 |
| 2005 | 8.284320 |
| 2006 | 7.873798 |
| 2007 | 8.375674 |
| 2008 | 8.874800 |
| 2009 | 7.572710 |
| 2010 | 9.043721 |
| 2011 | 8.610385 |
| 2013 | 9.013859 |
| 2015 | 7.623354 |
| 2017 | 7.852672 |
| 2018 | 7.363181 |
| 2019 | 8.357738 |
| 2020 | 8.830386 |
| 2021 | 8.642675 |
| 2022 | 7.975856 |
ggplot(standard_dev_all_839, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 839 average temperatures for water years 2005-2021
MK & SS for 839 (non-corrected)
sd_mk_839 <- mk.test(standard_dev_all_839$sd_2)
print(sd_mk_839)##
## Mann-Kendall trend test
##
## data: standard_dev_all_839$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -151.0000000 3461.6666667 -0.3247312sd_sens_839 <- sens.slope(standard_dev_all_839$sd_2)
print(sd_sens_839)##
## Sen's slope
##
## data: standard_dev_all_839$sd_2
## z = -2.5495, n = 31, p-value = 0.01079
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.051315272 -0.008246409
## sample estimates:
## Sen's slope
## -0.03122783Corrected
839 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_839_ad <- snotel_839_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_839_ad <- yearly_wy_aver_839_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_839_ad <- daily_wy_aver_839_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_839_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_839_ad <-daily_wy_aver_839_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_839_ad$date_temp_ad <- signif(daily_wy_aver2_839_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_839_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
839 SS (corrected)
standard_dev_839_ad <- daily_wy_aver_839_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_839_ad <- standard_dev_839_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_839_ad <- standard_dev_all_839_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_839_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 8.404486 |
| 1988 | 8.657820 |
| 1989 | 8.999163 |
| 1990 | 7.974877 |
| 1991 | 8.538551 |
| 1992 | 8.275130 |
| 1993 | 8.430766 |
| 1995 | 7.896968 |
| 1996 | 7.767745 |
| 1997 | 8.045426 |
| 1998 | 8.511594 |
| 1999 | 7.265049 |
| 2000 | 7.783278 |
| 2001 | 8.666438 |
| 2002 | 8.721243 |
| 2003 | 8.026530 |
| 2005 | 8.284037 |
| 2006 | 7.873931 |
| 2007 | 8.375662 |
| 2008 | 8.875360 |
| 2009 | 7.573403 |
| 2010 | 9.044246 |
| 2011 | 8.610915 |
| 2013 | 9.013498 |
| 2015 | 7.624141 |
| 2017 | 7.852753 |
| 2018 | 7.363589 |
| 2019 | 8.358178 |
| 2020 | 8.831049 |
| 2021 | 8.643116 |
| 2022 | 7.976653 |
ggplot(standard_dev_all_839_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 839 average temperatures for water years 1986-2021
MK & SS 839 (corrected)
sd_mk_839_ad <- mk.test(standard_dev_all_839_ad$sd_2)
print(sd_mk_839_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_839_ad$sd_2
## z = -0.16996, n = 31, p-value = 0.865
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -11.00000000 3461.66666667 -0.02365591sd_sens_839_ad <- sens.slope(standard_dev_all_839_ad$sd_2)
print(sd_sens_839_ad)##
## Sen's slope
##
## data: standard_dev_all_839_ad$sd_2
## z = -0.16996, n = 31, p-value = 0.865
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.02519708 0.02245395
## sample estimates:
## Sen's slope
## -0.002233571Upper San Juan 840
Morrisey 4/7/2004
snotel_840 <- SNOTEL_san_juan_Area %>%
filter(site_id == "840")#str(snotel_840) # check the date, usually a character.
snotel_840$Date <- as.Date(snotel_840$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_840_clean <- snotel_840 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_840_clean <- snotel_840_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_840_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_840_clean <- snotel_840_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_840_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_840_cull_count <- snotel_840_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_840_cull_count
# filtering for too few observations in a year
snotel_840_cull_count_days <- snotel_840_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_840_cull_count_days## # A tibble: 11 x 2
## # Groups: waterYear [11]
## waterYear n
## <dbl> <int>
## 1 1980 310
## 2 1981 274
## 3 1982 1
## 4 1983 318
## 5 1984 322
## 6 1985 275
## 7 1993 349
## 8 1994 345
## 9 2005 309
## 10 2006 325
## 11 2013 3141980, 1981, 1982, 1983, 1984, 1985, 1993, 1994, 2005, 2006, 2013 need to be culled.
snotel_840_clean_culled <- snotel_840_clean %>%
filter(waterYear >= "1986" & waterYear != "1993" & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_840_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_840_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_840_xts <- xts(snotel_840_clean_culled$temperature_mean, order.by = snotel_840_clean_culled$Date)
dygraph(temp_840_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_840_clean_culled <- snotel_840_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_840_xts <- xts(snotel_840_clean_culled$temperature_mean, order.by = snotel_840_clean_culled$Date)
dygraph(temp_840_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Upper San Juan 840
Morrisey 4/7/2004
snotel_840_adjusted <- snotel_840_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-04-07", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
840 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_840 <- snotel_840_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_840 <- yearly_wy_aver_840 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_840 <- daily_wy_aver_840 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_840$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_840 <-daily_wy_aver_840 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_840$date_temp <- signif(daily_wy_aver2_840$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_840, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
840 SD
standard_dev_840 <- daily_wy_aver_840 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_840 <- standard_dev_840 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_840 <- standard_dev_all_840 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_840 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1986 | 7.232619 |
| 1987 | 7.739156 |
| 1988 | 8.527104 |
| 1989 | 8.409877 |
| 1990 | 7.857130 |
| 1991 | 8.153466 |
| 1992 | 7.489499 |
| 1995 | 7.461187 |
| 1996 | 7.863786 |
| 1997 | 7.870552 |
| 1998 | 8.397044 |
| 1999 | 7.262936 |
| 2000 | 7.859027 |
| 2001 | 8.336002 |
| 2002 | 8.619541 |
| 2003 | 8.032941 |
| 2004 | 8.221773 |
| 2007 | 7.731123 |
| 2008 | 7.978427 |
| 2009 | 7.095790 |
| 2010 | 8.100698 |
| 2011 | 7.956902 |
| 2012 | 7.805592 |
| 2014 | 7.587053 |
| 2015 | 6.894972 |
| 2016 | 7.701017 |
| 2017 | 7.249233 |
| 2018 | 7.180902 |
| 2019 | 7.751539 |
| 2020 | 7.888215 |
| 2021 | 7.886812 |
| 2022 | 7.530559 |
ggplot(standard_dev_all_840, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 840 average temperatures for water years 2005-2021
MK & SS for 840 (non-corrected)
sd_mk_840 <- mk.test(standard_dev_all_840$sd_2)
print(sd_mk_840)##
## Mann-Kendall trend test
##
## data: standard_dev_all_840$sd_2
## z = -1.5406, n = 32, p-value = 0.1234
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -96.0000000 3802.6666667 -0.1935484sd_sens_840 <- sens.slope(standard_dev_all_840$sd_2)
print(sd_sens_840)##
## Sen's slope
##
## data: standard_dev_all_840$sd_2
## z = -1.5406, n = 32, p-value = 0.1234
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.031404215 0.001716915
## sample estimates:
## Sen's slope
## -0.0141416Corrected
840 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_840_ad <- snotel_840_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_840_ad <- yearly_wy_aver_840_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_840_ad <- daily_wy_aver_840_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_840_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_840_ad <-daily_wy_aver_840_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_840_ad$date_temp_ad <- signif(daily_wy_aver2_840_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_840_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
840 SS (corrected)
standard_dev_840_ad <- daily_wy_aver_840_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_840_ad <- standard_dev_840_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_840_ad <- standard_dev_all_840_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_840_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1986 | 6.729334 |
| 1987 | 7.211134 |
| 1988 | 7.975188 |
| 1989 | 7.874323 |
| 1990 | 7.334504 |
| 1991 | 7.648553 |
| 1992 | 6.993487 |
| 1995 | 6.945265 |
| 1996 | 7.326356 |
| 1997 | 7.346542 |
| 1998 | 7.819710 |
| 1999 | 6.756828 |
| 2000 | 7.298996 |
| 2001 | 7.766886 |
| 2002 | 8.021992 |
| 2003 | 7.460019 |
| 2004 | 7.506468 |
| 2007 | 7.729865 |
| 2008 | 7.977294 |
| 2009 | 7.094470 |
| 2010 | 8.099982 |
| 2011 | 7.955897 |
| 2012 | 7.804828 |
| 2014 | 7.585948 |
| 2015 | 6.894260 |
| 2016 | 7.700243 |
| 2017 | 7.247788 |
| 2018 | 7.180489 |
| 2019 | 7.751133 |
| 2020 | 7.887037 |
| 2021 | 7.885810 |
| 2022 | 7.530352 |
ggplot(standard_dev_all_840_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 840 average temperatures for water years 1986-2021
MK & SS 840 (corrected)
sd_mk_840_ad <- mk.test(standard_dev_all_840_ad$sd_2)
print(sd_mk_840_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_840_ad$sd_2
## z = 1.2811, n = 32, p-value = 0.2002
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 80.0000000 3802.6666667 0.1612903sd_sens_840_ad <- sens.slope(standard_dev_all_840_ad$sd_2)
print(sd_sens_840_ad)##
## Sen's slope
##
## data: standard_dev_all_840_ad$sd_2
## z = 1.2811, n = 32, p-value = 0.2002
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.005553565 0.027024754
## sample estimates:
## Sen's slope
## 0.0120862Vallecito 843 Morrisey 7/22/2005
snotel_843 <- SNOTEL_san_juan_Area %>%
filter(site_id == "843")#str(snotel_843) # check the date, usually a character.
snotel_843$Date <- as.Date(snotel_843$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_843_clean <- snotel_843 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_843_clean <- snotel_843_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_843_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_843_clean <- snotel_843_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_843_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_843_cull_count <- snotel_843_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_843_cull_count
# filtering for too few observations in a year
snotel_843_cull_count_days <- snotel_843_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_843_cull_count_days## # A tibble: 2 x 2
## # Groups: waterYear [2]
## waterYear n
## <dbl> <int>
## 1 1994 339
## 2 2003 3131994 & 2003 need to be culled.
snotel_843_clean_culled <- snotel_843_clean %>%
filter(waterYear != "1994" & waterYear != "2003")# & waterYear != "1994" & waterYear != "2005" & waterYear != "2006" & waterYear != "2013")# & waterYear != "2022") & waterYear != "2014") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_843_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_843_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_843_xts <- xts(snotel_843_clean_culled$temperature_mean, order.by = snotel_843_clean_culled$Date)
dygraph(temp_843_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_843_clean_culled <- snotel_843_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_843_xts <- xts(snotel_843_clean_culled$temperature_mean, order.by = snotel_843_clean_culled$Date)
dygraph(temp_843_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Vallecito 843
Morrisey 7/22/2005
snotel_843_adjusted <- snotel_843_clean_culled %>%
mutate(temp_ad = if_else(Date < "2005-07-22", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
843 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_843 <- snotel_843_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_843 <- yearly_wy_aver_843 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_843 <- daily_wy_aver_843 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_843$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_843 <-daily_wy_aver_843 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_843$date_temp <- signif(daily_wy_aver2_843$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_843, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
843 SD
standard_dev_843 <- daily_wy_aver_843 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_843 <- standard_dev_843 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_843 <- standard_dev_all_843 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_843 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.744307 |
| 1988 | 8.133975 |
| 1989 | 8.353802 |
| 1990 | 7.814664 |
| 1991 | 7.913458 |
| 1992 | 7.258470 |
| 1993 | 8.017543 |
| 1995 | 7.737939 |
| 1996 | 8.012757 |
| 1997 | 8.122829 |
| 1998 | 8.284299 |
| 1999 | 6.874497 |
| 2000 | 8.036954 |
| 2001 | 8.326535 |
| 2002 | 8.475624 |
| 2004 | 5.169906 |
| 2005 | 7.879202 |
| 2006 | 7.087987 |
| 2007 | 7.572592 |
| 2008 | 8.035842 |
| 2009 | 7.216959 |
| 2010 | 7.825626 |
| 2011 | 7.886538 |
| 2012 | 7.648815 |
| 2013 | 8.040098 |
| 2014 | 7.292347 |
| 2015 | 6.734206 |
| 2016 | 7.523001 |
| 2017 | 7.419976 |
| 2018 | 7.106676 |
| 2019 | 7.826098 |
| 2020 | 7.899858 |
| 2021 | 7.970432 |
| 2022 | 7.307690 |
ggplot(standard_dev_all_843, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 843 average temperatures for water years 2005-2021
MK & SS for 843 (non-corrected)
sd_mk_843 <- mk.test(standard_dev_all_843$sd_2)
print(sd_mk_843)##
## Mann-Kendall trend test
##
## data: standard_dev_all_843$sd_2
## z = -1.6307, n = 34, p-value = 0.103
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -111.000000 4550.333333 -0.197861sd_sens_843 <- sens.slope(standard_dev_all_843$sd_2)
print(sd_sens_843)##
## Sen's slope
##
## data: standard_dev_all_843$sd_2
## z = -1.6307, n = 34, p-value = 0.103
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.029639051 0.001693822
## sample estimates:
## Sen's slope
## -0.01383409Corrected
843 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_843_ad <- snotel_843_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_843_ad <- yearly_wy_aver_843_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_843_ad <- daily_wy_aver_843_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_843_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_843_ad <-daily_wy_aver_843_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_843_ad$date_temp_ad <- signif(daily_wy_aver2_843_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_843_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
843 SS (corrected)
standard_dev_843_ad <- daily_wy_aver_843_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_843_ad <- standard_dev_843_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_843_ad <- standard_dev_all_843_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_843_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.167709 |
| 1988 | 7.554905 |
| 1989 | 7.775142 |
| 1990 | 7.236968 |
| 1991 | 7.363376 |
| 1992 | 6.713681 |
| 1993 | 7.440691 |
| 1995 | 7.144630 |
| 1996 | 7.409097 |
| 1997 | 7.529924 |
| 1998 | 7.661187 |
| 1999 | 6.355448 |
| 2000 | 7.405969 |
| 2001 | 7.695693 |
| 2002 | 7.828232 |
| 2004 | 4.700776 |
| 2005 | 7.205286 |
| 2006 | 7.087285 |
| 2007 | 7.571307 |
| 2008 | 8.034938 |
| 2009 | 7.215649 |
| 2010 | 7.824301 |
| 2011 | 7.885342 |
| 2012 | 7.647841 |
| 2013 | 8.038906 |
| 2014 | 7.290741 |
| 2015 | 6.733530 |
| 2016 | 7.521405 |
| 2017 | 7.418681 |
| 2018 | 7.106001 |
| 2019 | 7.824924 |
| 2020 | 7.898803 |
| 2021 | 7.969378 |
| 2022 | 7.306728 |
ggplot(standard_dev_all_843_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 843 average temperatures for water years 1986-2021
MK & SS 843 (corrected)
sd_mk_843_ad <- mk.test(standard_dev_all_843_ad$sd_2)
print(sd_mk_843_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_843_ad$sd_2
## z = 1.601, n = 34, p-value = 0.1094
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 109.0000000 4550.3333333 0.1942959sd_sens_843_ad <- sens.slope(standard_dev_all_843_ad$sd_2)
print(sd_sens_843_ad)##
## Sen's slope
##
## data: standard_dev_all_843_ad$sd_2
## z = 1.601, n = 34, p-value = 0.1094
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.003122605 0.025938612
## sample estimates:
## Sen's slope
## 0.01173229Wolf Creek Summit 874
Morrisey 7/12/2004
snotel_874 <- SNOTEL_san_juan_Area %>%
filter(site_id == "874")#str(snotel_874) # check the date, usually a character.
snotel_874$Date <- as.Date(snotel_874$date) #change date from character to date format, capitalize to work with Water year functon from NWIS.
#THIS WILL CHANGE FOR EACH STATION
snotel_874_clean <- snotel_874 %>% # filter for the timeframe
filter(Date >= "1979-10-01" & Date <= "2022-09-30") %>%
#filter(temperature_mean >= -30 & temperature_mean <= 20) %>% # removing outliers
addWaterYear() %>%
mutate(daymonth = format(as.Date(Date), "%d-%m")) %>%
na.omit()
#adding water day using difftime (SUPER COOL. example from [this](https://stackoverflow.com/questions/48123049/create-day-index-based-on-water-year))
snotel_874_clean <- snotel_874_clean %>%
group_by(waterYear)%>%
mutate(waterDay = (as.integer(difftime(Date, ymd(paste0(waterYear - 1 ,'-09-30')), units = "days"))))# Check for outliers
ggplot(snotel_874_clean, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
snotel_874_clean <- snotel_874_clean %>%
mutate(temp_diff = abs(temperature_min - temperature_max)) %>%
filter(temperature_mean > -40) %>%
filter(temperature_mean < 25)ggplot(snotel_874_clean, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
Per Steven’s advice: :If there are more than 15 missing days, then…remove that year”
# filtering for temperature anomalies
#snotel_874_cull_count <- snotel_874_clean %>%
# filter(temperature_min > -40) %>%
# count(waterYear)
#snotel_874_cull_count
# filtering for too few observations in a year
snotel_874_cull_count_days <- snotel_874_clean %>%
group_by(waterYear) %>%
count(waterYear) %>%
filter(n < 350)
snotel_874_cull_count_days## # A tibble: 11 x 2
## # Groups: waterYear [11]
## waterYear n
## <dbl> <int>
## 1 1990 295
## 2 1994 313
## 3 2003 311
## 4 2006 340
## 5 2007 340
## 6 2008 315
## 7 2009 320
## 8 2010 303
## 9 2013 304
## 10 2016 338
## 11 2017 3141990, 1994, 2003, 2006, 2007, 2008, 2009, 2010, 2013, 2016, 2017 need to be culled.
snotel_874_clean_culled <- snotel_874_clean %>%
filter(waterYear != "1990" & waterYear != "1994" & waterYear != "2003" & waterYear != "2006" & waterYear != "2007" & waterYear != "2008" & waterYear != "2009" & waterYear != "2010" & waterYear != "2013" & waterYear != "2016" & waterYear != "2017") #%>%
#filter(temperature_mean > -49)
ggplot(snotel_874_clean_culled, aes(x = Date, y = temp_diff)) +
geom_point() + #lwd = 2) +
theme_few() +
#geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature varience (°C)') +
xlab('Date')
ggplot(snotel_874_clean_culled, aes(x = Date, y = temperature_mean)) +
geom_point() + #lwd = 2) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('Daily temperature (°C)') +
xlab('Date')
temp_874_xts <- xts(snotel_874_clean_culled$temperature_mean, order.by = snotel_874_clean_culled$Date)
dygraph(temp_874_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") #snotel_874_clean_culled <- snotel_874_clean_culled %>%
# filter(temperature_mean < 19.3)
temp_874_xts <- xts(snotel_874_clean_culled$temperature_mean, order.by = snotel_874_clean_culled$Date)
dygraph(temp_874_xts) %>%
dyAxis("y", label = "Daily mean temperature (°C)") Wolf Creek Summit 874
Morrisey 7/12/2004
snotel_874_adjusted <- snotel_874_clean_culled %>%
mutate(temp_ad = if_else(Date < "2004-07-12", ((5.3*10^(-7))*(temperature_mean^(4))+(3.72*10^(-5))*(temperature_mean^(3))-(2.16*10^(-3))*(temperature_mean^(2))-(7.32*10^(-2))*(temperature_mean)+1.37)+temperature_mean, temperature_mean))Non-corrected
874 Detrended
#using the clean culled df:
#average water year temperature
yearly_wy_aver_874 <- snotel_874_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp = mean(temperature_mean))
#Average temperature by day for all water years:
daily_wy_aver_874 <- yearly_wy_aver_874 %>%
group_by(daymonth) %>%
mutate(aver_day_temp = mean(aver_ann_temp))
#average mean temperature by day for the period of record:
daily_wy_aver_874 <- daily_wy_aver_874 %>%
group_by(daymonth) %>%
mutate(all_ave_temp = mean(daily_wy_aver_874$aver_day_temp))
# try to show all years as means.
daily_wy_aver2_874 <-daily_wy_aver_874 %>%
group_by(waterDay) %>%
mutate(date_temp = mean(temperature_mean))
daily_wy_aver2_874$date_temp <- signif(daily_wy_aver2_874$date_temp,3) #reduce the sig figs
ggplot(daily_wy_aver2_874, aes(x = waterDay, y = date_temp))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
874 SD
standard_dev_874 <- daily_wy_aver_874 %>%
group_by(waterYear) %>%
mutate(residual = (all_ave_temp-aver_ann_temp)+temperature_mean-aver_day_temp) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_874 <- standard_dev_874 %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_874 <- standard_dev_all_874 %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_874 %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.854854 |
| 1988 | 8.380441 |
| 1989 | 8.376871 |
| 1991 | 7.952272 |
| 1992 | 7.274045 |
| 1993 | 7.924599 |
| 1995 | 7.536845 |
| 1996 | 8.020204 |
| 1997 | 8.005351 |
| 1998 | 8.406529 |
| 1999 | 7.057906 |
| 2000 | 8.042917 |
| 2001 | 8.360075 |
| 2002 | 8.824249 |
| 2004 | 8.065839 |
| 2005 | 7.482529 |
| 2011 | 8.051681 |
| 2012 | 7.789860 |
| 2014 | 7.447867 |
| 2015 | 6.928269 |
| 2018 | 7.223986 |
| 2019 | 7.800112 |
| 2020 | 7.863897 |
| 2021 | 8.007803 |
| 2022 | 7.236522 |
ggplot(standard_dev_all_874, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Non-corrected standard deviation of SNOTEL 874 average temperatures for water years 2005-2021
MK & SS for 874 (non-corrected)
sd_mk_874 <- mk.test(standard_dev_all_874$sd_2)
print(sd_mk_874)##
## Mann-Kendall trend test
##
## data: standard_dev_all_874$sd_2
## z = -1.4247, n = 25, p-value = 0.1543
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## -62.0000000 1833.3333333 -0.2066667sd_sens_874 <- sens.slope(standard_dev_all_874$sd_2)
print(sd_sens_874)##
## Sen's slope
##
## data: standard_dev_all_874$sd_2
## z = -1.4247, n = 25, p-value = 0.1543
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.046813576 0.005791174
## sample estimates:
## Sen's slope
## -0.01681136Corrected
874 Detrended (corrected)
#using the clean culled df:
#average water year temperature
yearly_wy_aver_874_ad <- snotel_874_adjusted %>%
group_by(waterYear) %>%
mutate(aver_ann_temp_ad = mean(temp_ad))
#Average temperature by day for all water years:
daily_wy_aver_874_ad <- yearly_wy_aver_874_ad %>%
group_by(daymonth) %>%
mutate(aver_day_temp_ad = mean(aver_ann_temp_ad))
#average mean temperature by day for the period of record:
daily_wy_aver_874_ad <- daily_wy_aver_874_ad %>%
group_by(daymonth) %>%
mutate(all_ave_temp_ad = mean(daily_wy_aver_874_ad$aver_day_temp_ad))
# try to show all years as means.
daily_wy_aver2_874_ad <-daily_wy_aver_874_ad %>%
group_by(waterDay) %>%
mutate(date_temp_ad = mean(temp_ad))
daily_wy_aver2_874_ad$date_temp_ad <- signif(daily_wy_aver2_874_ad$date_temp_ad,3) #reduce the sig figs
ggplot(daily_wy_aver2_874_ad, aes(x = waterDay, y = date_temp_ad))+
geom_line(size= 0.7) +
theme_few() +
ylab('Average Daily temperature (°C)') +
xlab('Day of water year')
874 SS (corrected)
standard_dev_874_ad <- daily_wy_aver_874_ad %>%
group_by(waterYear) %>%
#filter(waterYear >= 1987 & waterYear <= 2021) %>%
mutate(residual = (all_ave_temp_ad-aver_ann_temp_ad)+temp_ad-aver_day_temp_ad) %>%
mutate(deviation = abs(residual-lag(residual)))
standard_dev_all_874_ad <- standard_dev_874_ad %>%
group_by(waterYear) %>%
mutate(nmbr = n())
standard_dev_all_874_ad <- standard_dev_all_874_ad %>%
group_by(waterYear) %>%
mutate(resid_mean = mean(residual)) %>%
mutate(sd_1 = residual-resid_mean) %>%
mutate(sd_2 = (((sum((sd_1)^2))/((nmbr-1))))^(0.5)) %>%
distinct(sd_2, .keep_all = TRUE) %>%
select(waterYear, sd_2)
standard_dev_all_874_ad %>%
kable(.,'html') %>%
kable_styling() %>%
scroll_box(width='250px',height='500px')| waterYear | sd_2 |
|---|---|
| 1987 | 7.327962 |
| 1988 | 7.851838 |
| 1989 | 7.848099 |
| 1991 | 7.458480 |
| 1992 | 6.775174 |
| 1993 | 7.395057 |
| 1995 | 7.010946 |
| 1996 | 7.470539 |
| 1997 | 7.481361 |
| 1998 | 7.826754 |
| 1999 | 6.568878 |
| 2000 | 7.465159 |
| 2001 | 7.779500 |
| 2002 | 8.206390 |
| 2004 | 7.444487 |
| 2005 | 7.482090 |
| 2011 | 8.052315 |
| 2012 | 7.790524 |
| 2014 | 7.447936 |
| 2015 | 6.928554 |
| 2018 | 7.224892 |
| 2019 | 7.800202 |
| 2020 | 7.864350 |
| 2021 | 8.008044 |
| 2022 | 7.236992 |
ggplot(standard_dev_all_874_ad, aes(x = waterYear, y = sd_2))+
geom_line(size= 0.7) +
theme_few() +
geom_smooth(method = "lm", se=FALSE) +
ylab('SD') +
xlab('Water year')
Morrisey-corrected standard deviation of SNOTEL 874 average temperatures for water years 1986-2021
MK & SS 874 (corrected)
sd_mk_874_ad <- mk.test(standard_dev_all_874_ad$sd_2)
print(sd_mk_874_ad)##
## Mann-Kendall trend test
##
## data: standard_dev_all_874_ad$sd_2
## z = 0.81742, n = 25, p-value = 0.4137
## alternative hypothesis: true S is not equal to 0
## sample estimates:
## S varS tau
## 36.000 1833.333 0.120sd_sens_874_ad <- sens.slope(standard_dev_all_874_ad$sd_2)
print(sd_sens_874_ad)##
## Sen's slope
##
## data: standard_dev_all_874_ad$sd_2
## z = 0.81742, n = 25, p-value = 0.4137
## alternative hypothesis: true z is not equal to 0
## 95 percent confidence interval:
## -0.006890724 0.033504267
## sample estimates:
## Sen's slope
## 0.006882736Figure captions are not correct.